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1.  BACKGROUND 


In  conjunction  with  a  study  of  the  benefits  of  insensitive  munition  (IM)  technology,  a  need  arose  to 
develop  a  methodology  for  prediction  of  the  probabilities  of  propagation  of  reaction  (detonation  or 
burning)  and  mechanical  damage  between  a  detonating  stack  of  ammunition  and  its  neighbors  as  functions 
of  the  distance  between  them.  Such  a  capability  could  provide  input  to  further  analyses  to  predict  losses 
in  a  variety  of  combat  scenarios. 

An  informal  survey  of  propagation  accident  reports  (summarized  in  Appendix  A)  indicates  that  the 
most  common  mechanism  of  reaction  propagation  involves  ignition  of  fires  (often  in  combustible 
packaging)  by  fragments,  debris,  or  firebrands  from  the  source  explosion  and  subsequent  violent  reaction 
of  munitions  in  those  fires.  The  resulting  chain  of  events  may  take  hours  or  even  days  to  unfold.  This 
scenario  appears  to  be  too  complex  and  variable  to  model  at  the  present  time. 

In  some  other  cases,  fragments  from  the  source  (donor)  explosion  can  promptly  damage,  ignite  mild 
to  violent  burning  in,  or  detonate  the  energetic  components  (not  including  any  combustible  packaging)  of 
munitions  in  nearby  (acceptor)  stacks.  When  detonation  results,  reaction  can  propagate  further  by  the 
same  mechanism,  rapidly  consuming  large  quantities  of  ammunitioa  When  only  burning  or  mechanical 
damage  results,  the  acceptor  stack  may  be  totally  or  partially  destroyed  but  no  further  propagation  by  the 
fragment  mechanism  ensues.  This  scenario  is  amenable  to  modeling. 

In  order  to  span  the  range  of  munition  vulnerability,  we  wanted  to  obtain  predictions  applicable  to 
typical  (thick-walled)  artillery  projectiles  and  (thin-walled)  missiles  stacked  on  pallets.  We  also  wanted 
to  obtain  predictions  for  a  single  artillery  projectile  donor  representing  an  attacking  munition  which  might 
be  required  to  start  the  donor-acceptor  chain  in  an  analysis.  We  have  chosen  palletized  and  single  M107 
155-mm  projectiles  and  palletized  TOW-2A  missiles  as  representative  items. 

2.  APPROACH 

We  determined  that  the  following  elements  are  required  to  successfully  model  propagation  among 
stacks  of  these  items: 

•  descriptions  of  the  stack  storage  arrangements 
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•  arena  test  data  for  donor  stack  munitions  describing  the  initial  fragment  mass,  velocity  and  shape 
distributions,  or  estimates  of  such  data 

•  a  treatment  of  fragment  trajectories  to  determine  hazard  probabilities 

•  descriptions  of  the  vulnerable  components  of  the  acceptor  stack  munitions 

•  criteria  for  the  initiation  of  detonation  and  burning  as  well  as  for  mechanical  damage. 

The  pallet  arrangements  commonly  used  for  the  weapons  of  interest  have  been  determined.  These 
form  the  building  blocks  for  larger  stacks. 

Arena  test  data  is  generally  available  only  for  warheads  whose  performance  is  measured  by 
fragmentation  statistics.  That  is,  data  is  available  for  the  M107  projectile  but  not  for  the  shaped-charge 
warhead  and  rocket  motors  found  in  the  TOW-2A.  We  had  to  develop  a  method  to  estimate  the  data  for 
the  latter  items. 

We  decided  to  modify  an  existing  computer  program  called  FRAGHAZ  (McClesky  1988)  in  order  to 
compute  fragment  trajectories  and  the  desired  probabilities.  FRAGHAZ  was  developed  to  predict  the 
hazard  to  a  human  target  due  to  fragmentation  from  an  exploding  ammunition  stack.  Fragmentation  data 
for  a  number  of  munitions  is  provided  with  the  program. 

Descriptions  of  the  vulnerable  components  of  these  weapons  were  sometimes  difficult  to  come  by. 
They  were  obtained  from  a  variety  of  sources.  In  some  cases,  best  guesses  were  used.  In  all  cases,  the 
geometries  were  simplified  to  simple  cylindrical  metal  shells  filled  with  energetic  materials. 

An  initiation  criterion  for  detonation  is  well  known  and  a  simple  criterion  for  burning  may  be 
developed  from  research  data  pertinent  to  this  phenomenon. 

3.  FRAGHAZ 

FRAGHAZ  uses  a  fourth  order  Runge-Kutta  scheme  to  compute  trajectories  for  each  fragment 
specified  in  the  arena  test  input  data.  The  trajectories  include  effects  of  ricochet  and  wind  (if  desired). 
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All  the  fragments  are  assumed  to  emanate  from  a  single  vertical  line  extending  from  the  base  to  the  top 
of  munitions  in  the  face  of  the  stack  as  shown  in  Figure  1.  This  defines  the  stack  axis.  The  initial 
fragment  height  is  selected  at  random  within  this  range.  The  initial  fragment  velocity  and  elevation  angle 
are  selected  at  random  in  a  range  near  the  values  given  in  the  input  data.  Parameters  used  to  determine 
the  fragment  drag  coefficient  are  determined  as  functions  of  a  randomly  selected  value.  The  program 
replicates  these  trajectory  computations  many  times  with  randomly  selected  environmental  conditions. 
While  FRAGHAZ  provides  "Monte  Carlo"  and  "Full  Factorial"  options,  the  random  selections  referred 
to  in  the  foregoing  description  apply  only  to  the  Monte  Carlo  option. 
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A  downrange  "hazard  volume,"  also  shown  in  Figure  1,  is  defined.  This  consists  of  a  cylindrical 
sector  originating  at  the  stack  axis  having  a  specified  angular  width  (equivalent  to  the  azimuthal  sector 
associated  with  arena  test  data  collection)  and  having  the  height  of  the  target  (nominally  a  5.72-ft  tall 
standing  man).  It  is  divided  into  a  number  of  100-ft  annular  segments.  As  each  fragment  passes  through 
a  hazard-volume  segment,  various  hit  probabilities  and  fragment  densities  are  accumulated  and  represented 
as  functions  of  the  downrange  distance  associated  with  the  midpoint  of  the  segment  A  determination  of 
whether  or  not  the  fragment  is  hazardous  based  on  its  kinetic  energy  is  made.  The  hazard  hit  probabilities 
depend  on  the  ratio  of  the  presented  area  of  the  target  to  the  presented  area  of  the  hazard-volume  segment 
with  respect  to  hazardous  fragments  as  shown  in  Figure  2. 


Figure  2.  FRAGHAZ  representation  of  target  (acceptor)  vulnerability  and  algorithm  for  determining 
hazard  probability  (from  NSWC  TR  87-59). 


4.  FRAGPROP 


While  its  basic  features  were  retained,  the  modifications  to  FRAGHAZ  were  extensive.  They 
proceeded  in  two  phases:  the  first  to  streamline  the  existing  code  and  the  second  to  implement  the 
required  new  models.  The  resulting  code  was  renamed  FRAGPROP.  A  listing  of  the  source  code  is 
provided  in  Appendix  B.  A  flowchart  applicable  to  both  FRAGHAZ  and  FRAGPROP  is  shown  in 
Figure  3. 
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The  original  program  made  minimal  use  of  modularization  and  was  unwieldy  to  work  with.  We  broke 
the  code  down  into  subroutines.  As  a  result,  the  overall  logic  is  easily  followed  in  the  main  program. 
We  changed  the  variable  names  to  more  descriptive  forms  which  conform  to  FORTRAN  naming 
conventions.  In  conjunction  with  this,  we  eliminated  most  of  the  type  declarations.  Many  of  the  variables 
and  computations  were  unnecessary  for  our  purposes  and  we  eliminated  them.  Some  expressions,  although 
correct,  were  recast  in  order  to  eliminate  mixed-type  operations.  Elseif  constructions  were  eliminated  for 
clarity  and  less  convoluted  logic  was  substituted  where  possible.  Execution  speed  was  increased  by 
moving  expressions  which  had  been  repeatedly  evaluated  within  loops  to  positions  outside  those  loops  and 
by  eliminating  redundant  computations.  For  simplicity,  we  also  eliminated  the  Full  Factorial  option. 


New  features  and  algorithms  were  also  incorporated.  Input  was  reorganized  to  read  separate  files  for 
the  run  parameters,  the  donor  description,  and  the  acceptor  description.  Fragmentation  data  input  was 
rewritten  to  reject  fragments  not  meeting  user-specified  criteria  for  minimum  mass  and  initial  elevation 
angle.  The  downrange  segment  size  (which  had  been  fixed  at  100  ft)  was  generalized  to  accept  user 
specificatioa  Computation  of  vulnerable  areas  for  die  acceptor  elements  and  algorithms  for  determination 
of  the  detonation,  bum,  and  mechanical  damage  lethality  of  each  fragment  including  the  effects  of 
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container  penetration  were  added.  Finally,  a  determination  and  report  of  the  minimum  range  at  which  the 
analysis  is  applicable  was  implemented.  These  algorithms  are  described  in  detail  in  the  following  sections. 

5.  DONOR  MODELS 


5.1  Stack  Descriptioa  The  donor  stack  is  described  as  in  the  original  FRAGHAZ  code.  The  height 
of  die  bottom  of  the  stack  and  the  height  of  munitions  in  the  stack  are  required  to  provide  bounds  for  the 
initial  height  of  each  fragment  The  size  of  the  stack  contributing  to  downrange  fragmentation  is  described 
by  specifying  the  number  of  "interaction  areas"  present  on  the  downrange  face  of  the  stack.  An  interaction 
area  is  formed  by  each  side-by-side  pair  of  munitions.  Thus,  an  arrangement  of  three  projectiles  in  a  row 
forms  two  interaction  areas.  If  two  such  arrangements  are  stacked  one  on  top  of  the  other,  four  interaction 
areas  result.  However,  if  two  such  arrangements  are  set  side  by  side  (six  rounds  in  a  row),  five  interaction 
areas  are  produced.  Interaction  areas  are  considered  instead  of  individual  munitions  because  the  effect  of 
focusing  along  the  plane  of  symmetry  between  the  rounds  augments  the  fragment  velocities.  Interaction 
areas  are  specified  along  with  arena  test  data  for  multiple  munitions.  However,  where  data  for  single 
munitions  (including  estimated  data)  must  be  used,  the  donor  unit  is  the  single  munition  and  the  number 
of  munitions  on  the  downrange  face  specifies  the  stack  size. 

5.2  Arena  Test  Fragmentation  Data.  Arena  test  data  forms  the  basis  for  describing  the  fragmentation 
characteristics  of  a  particular  munition.  Data  for  interacting  Ml 07  projectiles  including  parameters  for 
215  fragments  is  provided  with  FRAGHAZ.  Each  fragment  recovered  in  the  test  is  characterized  by  five 
parameters:  the  polar  zone  of  its  origin,  its  mass,  velocity,  and  area  to  mass  ratio  as  well  as  the  ratio  of 
its  maximum  to  average  presented  area  (to  account  for  the  effects  of  tumbling  on  aerodynamic  drag). 

Fragment  multipliers  (McClesky  1988)  are  used  to  normalize  the  data  for  scaling  to  any  stack  size 
(i.e.,  any  number  of  interaction  areas  or  munitions)  and  any  hazard  volume  azimuthal  sector  width.  One 
fragment  multiplier  must  be  supplied  for  each  polar  zone  represented  in  the  data.  The  formula  for 
computing  the  fragment  multiplier  for  the  rth  polar  zone,  4^,  is 


Mf^i) 


1 

NiVs(4>i)  ’ 
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where  Nj  is  the  number  of  interaction  areas  (or  single  munitions,  if  appropriate)  in  the  face  of  the  stack 
from  which  the  data  was  collected  and  Ys(<j>})  is  the  azimuthal  sector  width  over  which  fragments  were 
collected  in  the  ith  polar  zone. 

5.3  Adapting  Arena  Test  Fragmentation  Data.  Arena  test  data  for  single  munitions  tabulated  in  the 
Joint  Munitions  Effectiveness  Manual  of  Fragmentation  Data  (FM  101-62-3)  is  not  in  the  FRAGHAZ 
format  and  does  not  provide  all  of  the  parameters  required.  Data  for  individual  fragments  is  not  provided. 
Rather,  the  mean  fragment  weight  and  integrated  number  of  fragments  for  50-grain  weight  increments  are 
tabulated  along  with  the  initial  velocity  for  each  polar  zone.  A  shape  factor,  K  (=0.5126  cm^g2^  for  the 
Ml 07),  for  all  fragments,  is  given  and  can  be  used  in  the  expression, 


af 

avg 

"®T 


to  estimate  the  area  to  mass  ratio  of  a  fragment. 

We  expanded  this  data  to  generate  a  table  of  870  individual  fragments  (as  required  by  FRAGPROP) 
by  rounding  the  integrated  number  of  fragments  for  each  polar  zone  and  fragment  mass  increment  to  an 
integral  value  and  generating  table  entries  for  that  number  of  individual  fragments.  These  fragments  were 
each  assigned  the  mean  weight  of  the  associated  weight  increment  and  the  mean  velocity  of  the  associated 
polar  zone.  In  doing  this,  it  was  possible  to  eliminate  small  fragments.  Since  the  ratio  of  the  fragments’ 
maximum  to  average  presented  area  is  not  given,  we  used  the  average  value  taken  from  the  FRAGHAZ 
input  data  for  M107  stacks. 


af, 


avg 


1.479. 


For  this  arena  test,  there  is  only  one  unit  in  the  stack  and  the  integrated  number  of  fragments  is 
adjusted  to  apply  to  the  entire  360°  azimuthal  sector.  Thus,  the  fragment  multipliers  are  independent  of 
polar  zone  (MF  =  1/360  =  2.7778  x  10"3). 
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A  capability  to  invert  the  polar  distribution  of  the  donor  fragmentation  data  was  added  to  FRAGPROP 
in  order  to  allow  representation  of  a  nose-down  attacking  projectile  as  the  initial  event  in  a  propagation 
chain. 

5.4  Estimating  Fragmentation  Data.  The  TOW-2A  missile  is  not  a  fragmenting  weapon  and 
fragmentation  data  is  not  available.  We  needed  to  estimate  data  having  a  representative  distribution  of 
values  of  the  five  parameters  required  by  the  FRAGPROP  input  Several  analytical  techniques  applicable 
to  single  cylindrical  munitions  are  available  for  this  purpose.  The  fragmentation  unit  in  this  case  is  a 
single  munition  rather  than  an  interaction  area.  Thus,  focusing  effects  are  ignored  and  the  number  of 
rounds  is  specified  in  lieu  of  the  number  of  interaction  areas.  Geometric  data  required  to  represent  the 
TOW-2A  warhead  and  rocket  motor  as  cylindrical  components  and  energetic  material  performance  data 
as  needed  for  use  with  the  following  analyses  are  summarized  in  Appendix  C. 

The  Mott  equation  (see  Victor  1994)  describes  the  distribution  of  fragment  mass.  The  fraction,  f,  of 
fragments  having  a  mass  greater  than  mf  is  given  by 


f(mf)  =  ^ex p(-2mf/m^)  , 


where  mf  is  the  average  fragment  mass.  This  equation  may  be  solved  for  mf. 


mf(f)  =  -i.m71n(f2)  . 


By  selecting  random  values  of  f  which  are  evenly  distributed  between  0  and  1,  this  equation  may  be  used 
to  generate  sets  of  fragment  masses  having  the  Mott  distribution.  The  total  number  of  fragments  is  given 
by 


nf  = 
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where  mc,  the  total  mass  of  all  fragments,  may  be  identified  with  the  casing  mass.  Victor  also  gives  the 
average  fragment  mass  as 


mf  =  —  \j  1  +p/2 
Pd 


where  A  is  a  constant,  pd  is  the  detonation  pressure,  p  is  the  ratio  of  the  mass  of  the  casing  to  the  mass 
of  the  charge,  hc  is  the  thickness  of  the  casing,  and  dx  is  the  diameter  of  the  explosive  or  propellant  charge 
(equal  to  the  inside  diameter  of  the  casing).  A  =  676.2  g-kbar/in3/2  for  mass  in  g,  pressure  in  kbar  and 
distance  in  inches.  The  Kamlet-Jacobs  formula  may  be  used  to  estimate  the  detonation  pressure  (in  GPa) 
if  it  is  not  otherwise  known. 


pd  =  l-SSSp^N^MQ-  , 


where  p0  is  the  unreacted  density  in  g/cm3,  N  is  the  number  of  moles  of  detonation  product  per  gram  of 
unreacted  material,  M  is  the  average  molecular  weight  of  the  detonation  product  gases,  and  Q  is  the 
chemical  energy  of  the  detonation  reaction  in  cal/g. 

No  method  exists  for  predicting  a  distribution  of  fragment  velocities.  However,  the  Gumey  equation 
(see  Dehn  1984)  may  be  used  to  estimate  a  single  velocity  for  all  fragments. 

s/2E 

vf  =  v  , 

vfa  +  1/2 


where  v/2E  is  a  property  of  the  energetic  material  (having  the  units  of  velocity)  known  as  the  Gumey 
constant.  Its  value  may  be  found  in  the  literature  or  estimated  (in  milliseconds)  using  the  equation  of 
Kamlet  and  Finger. 

,/2E  -  233^.  , 
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for  density  in  g/cm3  and  detonation  pressure  in  kbar.  Although  this  approach  gives  only  one  value  of  the 
fragment  velocity,  randomization  over  a  narrow  range  takes  place  within  FRAGPROP. 

For  estimating  the  area  to  mass  ratio  of  a  fragment,  Victor  (1994)  gives  the  shape  factor  as 
K  =  0.5199  cn^/g2^3  for  randomly  shaped  fragments.  We  are  not  aware  of  any  method  to  estimate  die 
ratio  of  a  fragment’s  maximum  to  average  presented  area,  so  we  assumed  a  fixed  value  of  1.5. 

The  TOW-2A  storage  orientation  is  horizontal  in  contrast  to  that  of  the  M107.  Thus,  it  is  necessary 
to  interpret  the  (rather  narrow)  polar  distribution  of  fragments  with  respect  to  the  munition  axis  as  an 
azimuthal  distribution  with  respect  to  the  stack  axis.  This  is  assumed  to  produce  an  azimuthal  sector  width 
of  only  20°.  Because  the  missile  components  are  treated  as  horizontal  cylinders,  the  polar  distribution 
with  respect  to  the  stack  axis  is  assumed  uniform.  Since  the  azimuthal  sector  width  is  independent  of 
polar  zone,  the  fragment  multiplier  is  a  constant  (MF  =  1/20  =  0.05). 

We  wrote  a  short  program  (incorporating  the  FRAGHAZ  random  number  generator)  to  produce 
FRAGPROP  input  data  for  a  general  warhead  and  rocket  motor  combination  in  accordance  with  the 
foregoing  algorithms  and  used  it  to  generate  data  for  365  fragments  from  the  TOW-2A  warhead  and  flight 
motor.  These  components  produce  both  aluminum  and  steel  fragments.  We  neglected  fragments  from 
the  launch  motor,  which  has  a  small  diameter.  The  program  is  listed  in  Appendix  D. 

6.  ACCEPTOR  MODELS 

6.1  Stack  Description.  The  basic  acceptor  stack  description  is  given  by  providing  dimensions  of 
height  (HJ,  width  (Wa),  and  depth  (Da).  The  storage  orientation  may  be  vertical,  as  with  the  M107,  or 
horizontal,  as  with  the  TOW-2A.  The  vulnerable  components  of  an  acceptor  may  include  a  warhead, 
rocket  motor,  or  both. 

6.2  Vulnerable  Areas.  Vulnerable  areas  at  the  front  and  top  of  the  stack  are  required.  Side  and  back 
areas  are  excluded  as  the  acceptor  stack  is  assumed  to  face  the  donor.  These  depend  on  the  weapon 
dimensions  and  the  storage  orientation  (horizontal  or  vertical). 

For  die  purpose  of  determining  hit  probabilities  the  total  front  and  top  areas  of  the  stack  are  considered 
vulnerable.  The  front  vulnerable  area  is 
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and  the  top  vulnerable  area  is 


Afh  -  HaWa 


Ath  =  DaWa  • 

For  mechanical  damage,  the  vulnerable  areas  are  the  areas  presented  to  the  front  and  top  faces  of  the 
stack  of  the  entire  weapon  (represented  as  a  cylinder  of  diameter  Dw  and  length  L^,)  multiplied  by  the 
number  of  weapons  in  that  face  of  the  stack  (Nf  or  Nt).  The  vulnerable  areas  depend  on  the  storage 
orientation.  For  vertical  storage,  the  front  vulnerable  area  is 

Afm  =  NfDwLw  * 

and  the  top  vulnerable  area  is 

\m  *  ».*  K«f  ■ 


For  horizontal  storage,  the  front  vulnerable  area  is  unchanged. 


Afm  “  NfDwLw  > 

while  the  top  vulnerable  area  is 

Atm  =  NtDwLw  . 


Similarly,  for  detonation  and  burning,  the  maximum  vulnerable  areas  are  those  of  the  energetic 
materials  in  the  warhead  and/or  rocket  motor  (cylinders  of  diameter  Dx  and  length  Lx)  presented  to  the 
front  and  top  faces  of  the  stack.  For  vertical  storage,  the  front  maximum  vulnerable  area  is 


Afx„„  =  NfDXLX 


and  the  top  maximum  vulnerable  area  is 


tx_ 
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For  horizontal  storage,  the  front  maximum  vulnerable  area  is 


Afx 


max 


=  NfDxLx  , 


and  the  top  maximum  vulnerable  area  is 


tx„ 


—  NtDxLx  . 


Here,  the  subscript,  x,  represents  either  the  warhead  explosive  or  rocket  motor  propellant 


With  respect  to  a  particular  fragment,  the  vulnerable  areas  for  detonation  and  burning  may  be  further 
limited  depending  on  the  velocity  and  diameter  of  that  fragment  Below  some  critical  threshold,  none  of 
the  presented  area  of  an  energetic  component  is  vulnerable.  Above  the  threshold,  part  or  all  of  the  areas 
is  vulnerable.  Partial  vulnerability  is  generally  a  function  of  the  maximum  obliquity,  0raax,  from  the 
normal  to  the  surface  of  the  weapon  component  at  which  an  impacting  fragment  can  produce  reaction. 
Expressions  for  0max  for  detonation  and  burning  are  developed  in  the  following  sections. 


6.3  Detonation  Vulnerability.  The  vulnerability  of  a  weapon  component  to  initiation  of  detonation 
by  fragment  impact  is  described  by  the  Jacobs-Roslund  formula  (Liddiard  and  Roslund  1993)  for  critical 
impact  velocity.  This  formula  applies  to  cylindrical  projectiles  having  diameters  greater  than  the  failure 
diameter  of  the  energetic  material. 


where  ajj.,  b^,  and  Cjr  are  characteristic  constants  which  have  been  determined  for  a  number  of  explosives. 
Values  of  the  Jacobs-Roslund  constants  used  in  conjunction  with  our  analysis  are  given  in  Appendix  C. 
It  is  assumed  that  the  fragment  always  strikes  the  vulnerable  component  with  an  average  value  of 
presented  area.  Thus,  the  fragment  diameter,  df,  may  be  determined  from  the  average  fragment  presented 
areas  assuming  a  circular  cross-section, 
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df  =  2 


avg 


M  n 


The  angle,  0,  is  the  obliquity  of  the  fragment  with  respect  to  the  normal  to  the  surface  of  the  casing.  The 
obliquity  dependence  of  the  critical  velocity  may  be  expressed  as 


v0) 


Vjr(Q) 
cos  9 


We  define  the  maximum  obliquity  that  can  produce  detonation  for  a  fragment  with  velocity,  vf,  such  that 

V0max)  =  vf  • 


Thus, 


cos0, 


max 


Vjr(O) 

vf 


and  0,,^  is  defined  for  vf  >  Vjr(0).  This  represents  the  critical  condition.  This  model  has  been  calibrated 
for  steel  fragments.  However,  we  also  used  this  calibration  for  the  aluminum  fragments  produced  by  the 
TOW-2A  warhead  casing.  As  a  worst  case,  we  assumed  that  the  fragments  were  flat-ended  cylinders 
(bjr  =  0  with  an  appropriate  value  for  Cjr)  and  that  the  critical  diameter  of  the  energetic  material  was 
always  smaller  than  the  diameter  of  the  fragment  Since  actual  fragments  rarely  exhibit  this  morphology, 
a  more  realistic  simulation  might  result  from  a  random  selection  of  the  critical  velocity  between  limits  for 
flat-  and  round-ended  cylindrical  fragments. 

6.4  Burning  Vulnerability.  The  vulnerability  of  a  weapon  component  to  initiation  of  burning  may 
be  associated  with  perforation  of  the  casing  by  the  fragment  This  is  particularly  hue  in  the  case  of  thin- 
walled  munitions  (Gilman).  Some  critical  residual  velocity,  vre,  may  be  required  to  produce  reaction.  The 
THOR  velocity  equation  relates  the  residual  velocity,  vr  to  the  initial  velocity,  vf,  for  a  penetrating 
fragment. 
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vr  =  vf  -10ay  (hcafavg)bvmfCvvfdv/(cos0)ev  . 


Again,  it  is  assumed  that  the  fragment  strikes  the  vulnerable  component  with  an  average  value  of  presented 
area.  With  vr  =  vrc  at  0  =  0,^, 


COS0, 


max 


10 


iav 


(Mu)b 


|l/ev 


mf  \f 


Vf-V; 


rc 


and  0max  is  defined  for 


vf  -vrc  >10av(hcaftvJbv  mfCv  vfdv . 

This  expression  represents  the  critical  condition.  Values  of  the  THOR  constants  used  in  conjunction  with 
our  analysis  are  given  in  Appendix  C.  Since  appropriate  values  of  the  critical  residual  velocity  are  not 
known,  we  used  vrc  =  0  as  a  worst  case.  This  model  has  also  been  calibrated  for  steel  fragments  but, 
again,  we  used  the  same  calibration  for  the  aluminum  fragments  produced  by  the  TOW-2A  warhead 
casing. 

6.5  Comparison  of  Detonation  and  Burning  Thresholds.  It  is  tempting  to  presume  that  the  less  violent 
burning  response  is  easier  to  produce  than  detonation.  However,  comparison  of  the  critical  velocities 
predicted  by  these  models  shows  that  this  is  not  always  the  case.  Critical  velocities  for  detonation  and 
burning  of  a  typical  fragment  are  plotted  as  functions  of  casing  thickness  in  Figure  4.  For  thin,  easily 
perforated  casings  (e.g.,  -1  mm  for  the  TOW-2A),  burning  is  produced  by  fragments  having  much  lower 
velocities  than  those  required  to  produce  detonation.  (This  result  would  still  hold  if  a  nonzero  critical 
residual  velocity  of  expected  magnitude  were  used.)  For  thicker  casings  (e.g.,  -15  mm  for  the  M107), 
casing  perforation  is  extremely  difficult  to  achieve  but  detonations  may  be  produced  at  lower  velocities 
due  to  the  transmitted  shock. 
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CASING  THICKNESS  (mm) 

Figure  4.  Comparison  of  critical  velocities  for  detonation  and  burning. 

The  point  of  intersection  of  the  two  threshold  curves  varies  with  the  shape  of  the  fragment 
(represented  by  the  area  to  mass  ratio).  For  ratios  lower  than  the  10  in2/lb  value  used  in  Figure  4  the 
point  of  intersection  is  shifted  toward  high  values  of  casing  thickness.  Thus,  results  are  sensitive  to 
assumptions  regarding  fragment  orientation  (presented  area)  at  impact. 

6.6  Vulnerable  Area  Reduction.  If  the  angle  subtended  at  the  center  of  the  cylinder  representing  a 
vulnerable  component  by  the  region  on  the  surface  over  which  an  impacting  fragment  is  lethal  with  respect 
to  detonation  or  burning  (0  <  0max)  is  2|3max  as  shown  in  Figure  5,  then  the  reduction  in  the  maximum 
vulnerable  area  of  that  component  is  given  by 

A  -  sinPmaxAmax  • 
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CYLINDRICAL 

VULNERABLE 

COMPONENT 


The  relationship  between  0max  and  pmax  depends  on  the  storage  orientation.  For  horizontal  storage, 

Pmax  =  ®max 


as  shown  in  Figure  6.  For  vertical  storage,  the  elevation  angle  of  the  fragment  trajectory,  otj,  must  be 
accounted  for,  as  shown  in  Figure  7.  Thus, 


cos  ft 


max 


cos0„ 


COS  Of 


6.7  Mechanical  Damage  Vulnerability.  A  single  mechanical  damage  criterion  based  on  kinetic  energy 
applies  to  all  weapon  components.  The  stack  is  considered  vulnerable  to  damage  when  the  kinetic  energy 
of  a  fragment  exceeds  a  specified  value,  regardless  of  the  obliquity.  Accurate  threshold  values  are  not 
known,  so  we  arbitrarily  chose  400  and  50  ft-lb  as  the  critical  kinetic  energy  for  the  Ml  07  and  TOW-2A, 
respectively.  The  smaller  value  for  the  TOW-2A  reflects  its  much  thinner  skin. 
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een  the  angular  region  of 


fragment  elevation  angle  for  vertical  storage. 
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7.  CONTAINER  PENETRATION 


The  TOW-2A  is  packaged  in  a  steel  container,  and  provisions  were  made  in  FRAGPROP  for 
estimating  the  reduction  in  fragment  mass  and  velocity  associated  with  penetration  of  the  container.  This 
occurs  both  when  fragments  are  launched  from  a  donor  stack  and  when  they  impact  an  acceptor  stack. 
In  addition  to  the  THOR  velocity  equation,  the  THOR  mass  equation  is  required. 

mr  =  mf-10am|hcaf^JbmrnfCm  vfdm/(cos0)e,n  . 

Use  of  the  THOR  equations  with  0  =  0  yields  the  maximum  residual  velocity  and  mass,  representing  the 
worst-case  scenario.  The  residual  diameter  is  determined  from  the  residual  mass  assuming  that  the  area 
to  mass  ratio  remains  unchanged. 

8.  MINIMUM  RANGE 

The  FRAGHAZ  analysis  is  applicable  as  long  as  the  angle  subtended  by  the  acceptor  stack  does  not 
exceed  the  width,  \|/,  of  the  azimuthal  sector  associated  with  the  donor  fragmentation  data.  The  minimum 
range  is  thus  given  by 


W„ 


*nun 


2tan(\[t/2) 


Note  that  rmin  is  zero  for  sectors  wider  than  180°. 

9.  COMPUTATIONAL  CONFIGURATIONS  AND  RESULTS 

9.1  Theater  of  Operations  Considerations.  Storage  regulations  applicable  to  basic  load  ammunition 
holding  areas  in  theaters  of  operations  limit  the  explosive  quantity  in  any  stack  to  4,000  kg  (8,818  lb). 
Such  stacks  must  be  separated  by  at  least  77  m  (253  ft).  (See  Army  Regulation  385-64.)  We  considered 
M107  and  TOW-2  A  stacks  containing  approximately  this  maximum  weight. 
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Each  M107  projectile  contains  15.4  lb  of  Composition  B.  Thus,  a  pallet  of  eight  projectiles  contains 

123.2  lb  and  a  stack  of  72  pallets  contains  8,870.4  lb. 


A  TOW-2A  missile  contains  6.7  lb  of  LX-14  plus  0.2  lb  of  other  high  explosives  in  the  warhead, 

1.2  lb  of  M7  propellant  in  the  launch  motor  and  7.0  lb  of  GCV  propellant  in  the  flight  motor  for  a  total 
of  15.1  lb.  A  pallet  of  12  missiles  contains  181.2  lb,  and  a  stack  of  48  pallets  contains  8,697.6  lb. 

9.2  Stack  Arrangements.  The  arrangement  of  pallets  affects  the  lethality  and  vulnerability  of  a  stack. 
Consideration  of  the  least  and  most  lethal  arrangements  as  well  as  the  least  and  most  vulnerable 
arrangements  provides  scope  to  a  determination  of  propagation  probabilities.  Pertinent  stack  parameters 
may  be  estimated  from  pallet  dimensions  given  in  Appendix  C. 

Each  M107  pallet  contains  eight  vertical  projectiles  in  an  arrangement  that  is  one  projectile  high,  four 
projectiles  wide,  and  two  projectiles  deep  (1  x  4  x  2).  The  72  pallets  may  be  arranged  in  any  permutation 
of  3  x  4  x  6,  where  the  first  dimension  is  the  height  of  the  stack  in  pallets,  the  second  dimension  is  the 
width  of  the  stack  in  pallets,  and  the  third  dimension  is  the  depth  of  the  stack  in  pallets.  Individual  pallets 
retain  their  1x4x2  orientations  for  all  permutations.  Table  1  gives  the  number  of  interaction  areas, 
stack  dimensions,  and  total  front  and  top  vulnerable  area  associated  with  each  permutation. 

The  4x3x6  pallet  arrangement  is  both  least  lethal  (having  the  fewest  interaction  areas)  and  least 
vulnerable  (having  the  smallest  vulnerable  area,  while  the  4  x  6  x  3  arrangement  is  both  most  lethal  and 
most  vulnerable. 


Each  TOW-2A  pallet  contains  12  horizontal  missiles  in  an  arrangement  that  is  3  missiles  high, 
1  missile  wide,  and  4  missiles  deep  (3  x  1  x  4).  The  48  pallets  may  be  arranged  in  any  permutation  of 
3x4x4.  Table  2  gives  the  number  of  missiles  on  the  front  face,  stack  dimensions,  and  total  front  and 
top  vulnerable  areas  associated  with  each  permutation. 

Based  on  the  number  of  missiles  on  the  front  face,  the  4  x  4  x  3  pallet  arrangement  is  most  lethal. 
The  other  two  arrangements  are  nominally  equal  in  lethality.  The  3x4x4  arrangement  is  most 
vulnerable,  while  the  4  x  3  x  4  arrangement  is  least  vulnerable.  The  missile  stacks  are  considerably  larger 
and  more  vulnerable  to  being  hit  than  the  artillery  projectile  stacks. 
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Table  1.  Lethality  and  Vulnerability  Parameters  for  M107  Stacks 


Pallet 

Arrangement 
(HX^x  nd) 

Number  of 
Interaction 
Areas  on 
Face  of  Stack 

Stack 

Dimensions 

Vulnerable 

Areas 

(front  +  top) 
(ft2) 

■el 

width 

(ft) 

depth 

(ft) 

3x4x6 

45 

8.00 

9.63 

7.20 

146.3 

3x6x4 

69 

8.00 

14.50 

4.77 

185.2 

4x3x6 

44 

10.67 

7.20 

7.20 

128.6 

4x6x3 

92 

10.67 

14.50 

3.55 

206.2 

6x3x4 

66 

16.00 

7.20 

4.77 

149.5 

6x4x3 

90 

16.00 

9.63 

3.55 

188.3 

Table  2.  Lethality  and  Vulnerability  Parameters  for  TOW-2A  Stacks 


Pallet 

Arrangement 
(rib  x  r^  x  i^) 

Number  of 
Missiles  on 
Face 
of  Stack 

Stack 

Dimensions 

Vulnerable 

Areas 

(front  +  top) 
(ft2) 

height 

(ft) 

width 

(ft) 

mm 

3x4x4 

36 

iiiia 

20.25 

14.64 

500.0 

4x3x4 

36 

15.17 

14.64 

425.3 

4x4x3 

48 

20.25 

10.96 

493.4 

The  single  M107  projectile  donor  may  be  erect  or  inverted  and  may  detonate  at  any  specified  height 
above  the  ground.  The  arrangements  producing  the  least  and  greatest  lethality  are  not  apparent.  However, 
preliminary  computations  indicated  that  the  inverted  round  at  the  surface  may  be  considered  least  lethal 
and  die  erect  round  at  a  50-ft  elevation  most  lethal. 

9.3  Pmnagatinn  Probabilities  and  Distances.  The  graphics  program  developed  in  conjunction  with 
FRAGPROP  plots  the  probabilities  of  exceeding  the  criteria  for  propagation  of  detonation  (labeled  D), 
burning  (B),  and  mechanical  damage  (M),  as  well  as  the  overall  (lethal  and  nonlethal)  hit  probability  (H) 
on  a  logarithmic  scale  as  a  function  of  range.  A  typical  example  is  shown  in  Figure  8.  The  dashed 
vertical  line  near  the  probability  axis  indicates  the  minimum  range  for  which  the  analysis  is  accurate. 
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Ranqe  (ft) 


Figure  8.  Probabilities  of  detonation,  burning,  mechanical  damage,  and  hit  as  functions  of  ranee  for  an 
M107  donor  stack  against  a  TOW-2A  acceptor  stack. 

Generally,  the  probabilities  decrease  with  increasing  range  and,  theoretically,  they  can  be  expected  to 
vanish  entirely  beyond  some  maximum  range.  For  comparison  purposes,  it  is  often  useful  to  consider  the 
range  at  which  a  probability  drops  to  some  small  value  (e.g„  1%).  The  small  probabilities  predicted  at 
long  range  may  become  nonmonotonic  prior  to  vanishing.  The  ability  to  predict  nonzero  probabilities  at 
long  range  is  limited  by  the  number  of  fragments  in  the  donor  data  and  the  number  of  replications. 
Generally,  the  results  are  not  accurate  near  vanishing  probability.  The  number  of  replications  required  to 
produce  monotonically  decreasing  probabilities  varies  inversely  with  the  number  of  fragments  describing 
the  donor. 
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We  made  FRAGPROP  computations  using  zero-wind  conditions  and  soil  constants  varying  between 
0.5  and  4.0  for  each  of  the  three  donors  (single  M107  projectiles,  stacked  M107  pallets,  and  stacked 
TOW-2A  pallets)  vs.  both  of  the  acceptors  (stacked  M107  and  TOW-2A  pallets).  Preliminary 
computations  were  made  to  determine  the  maximum  range  required  in  each  problem.  Input  conditions 
for  the  final  computations  are  summarized  in  Table  3.  These  values  depend  on  the  donor  and  apply  to 
both  acceptors. 


Table  3.  FRAGPROP  Computation  Input  Conditions 


Donor 

Number  of 
Fragments 

Number  of 
Replications 

Segment 

Size 

(ft) 

Maximum 

Range 

(ft) 

Single  M107 

870 

50 

50 

2,400 

M107  Stack 

215 

200 

100 

3,600 

TOW-2A  Stack 

365 

125 

20 

720 

In  order  to  determine  the  scope  of  probability  values,  two  computations  were  made  for  each  donor- 
acceptor  pair:  one  for  the  least  lethal  donor  stack  arrangement  against  the  least  vulnerable  acceptor  stack 
arrangement  and  one  for  the  most  lethal  donor  stack  arrangement  against  the  most  vulnerable  acceptor 
stack  arrangement  The  graphics  capability  allows  the  four  probabilities  to  be  plotted  as  functions  of  range 
for  two  problems  at  a  time  as  shown  in  Figures  9-14.  Associated  pairs  of  probability  curves  are  joined 
by  shading  lines. 

Results  for  single  M107  donors  against  palletized  M107  acceptors  are  shown  in  Figure  9.  The 
probabilities  decrease  gradually  with  range.  The  distance  at  which  the  detonation  propagation  probability 
drops  to  1%  (on  the  upper  curve  of  the  pair)  is  125  ft,  while  the  1%  distance  for  burning  is  only  1 13  ft. 
Thus,  burning  propagation  appears  unlikely  in  this  configuration.  The  probability  of  mechanical  damage 
remains  above  1%  to  a  range  of  687  ft,  while  the  hit  probability  drops  to  1%  at  1,369  ft. 
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Donor:  INVERTED  M107  0  0  ft  Donor:  ERECT  M107  @  50  ft 

Receptor  :  -4x3x6  M107  PALLETS  Acceptor  :  -4x6x3  M107  PALLETS 


Donor:  INVERTED  M107  e  0  ft  Donor:  ERECT  M107  @  50  ft 

Acceptor :  -4x3x4  T0W-2A  PALLETS  Acceptor :  3x4x4  T0W“2A  PALLETS 


donor  projectiles  against  TOW-2A  acceptor  stacks. 
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Donor:  4x3x6  M107  PALLETS  Donor: 4x6x3  M107  PALLETS 

Acceptor:  4x3x6  M107  PALLETS  Acceptor :  4x6x3  M107  PALLETS 


Figure  11.  Probabilities  of  detonation,  burning,  mechanical  damage,  and  hit  as  functions  of  range  for 
M107  donor  stacks  against  M107  acceptor  stacks. 


Donor:  4x3x6  M107  PALLETS  Donor:  4x6x3  M107  PALLETS 

Acceptor:  4x3x4  T0W-2A  PALLETS  Acceptor : 3x4x4  T0W-2A  PALLETS 


Range  (ft) 

Figure  12.  Probabilities  of  detonation,  burning,  mechanical  damage,  and  hit  as  functions  of  range  for 
M107  donor  stacks  against  TOW-2A  acceptor  stacks. 
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Donor:  3x4x4  T0W-2A  PALLETS  Donor:4x4x3  T0W-2A  PALLETS 

Acceptor:  4x3x6  M107  PALLETS  Acceptor :  4x6x3  M107  PALLETS 


Ranqe  (ft) 

Figure  13.  Probabilities  of  mechanical  damage  and  hit  as  functions  of  range  for  TOW-2A  donor  stacks 
against  Ml 07  acceptor  stacks. 


Donor:  3x4x4  T0W-2A  PALLETS  Donor: 4x4x3  T0W-2A  PALLETS 

Acceptor :  4x3x4  T0W-2A  PALLETS  Acceptor  :  3x4x4  T0W-2A  PALLETS 


Ranqe  (ft) 

Figure  14.  Probabilities  of  burning,  mechanical  damage,  and  hit  as  functions  of  range  forTOW-2A  donor 
stacks  against  TOW-2A  acceptor  stacks. 
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Similar  results  for  single  M107  donors  against  palletized  TOW-2A  acceptors  are  shown  in  Figure  10. 
The  1%  distance  for  detonation  is  162  ft,  while  a  bum  probability  greater  than  1%  persists  to  a  distance 
of  471  ft  The  probability  of  mechanical  damage  drops  to  1%  at  666  ft  The  range  over  which  the  overall 
hit  probability  remains  above  1%  is  1,551  ft 

Results  for  palletized  M107  donors  against  palletized  M107  acceptors  are  shown  in  Figure  11. 
Probabilities  near  unity  persist  to  a  greater  range.  This  arrangement  produces  a  1%  probability  of 
detonation  at  642  ft  while  burning  reactions  propagate  with  a  probability  of  1%  or  greater  over  distances 
to  685  ft  The  probability  of  mechanical  damage  remains  above  1%  to  2,871  ft.  The  1%  distance  for  the 
overall  hit  probability  is  2,950  ft. 

Similar  results  for  palletized  M107  donors  against  palletized  TOW-2A  acceptors  are  shown  in 
Figure  12.  This  arrangement  produces  a  1%  probability  of  detonation  at  643  ft,  while  burning  reactions 
propagate  with  a  probability  of  1%  or  greater  over  distances  to  1,715  ft.  A  probability  of  mechanical 
damage  greater  than  1%  persists  to  2,779  ft.  The  1%  distance  for  the  overall  hit  probability  is  3,080  ft 

Results  for  palletized  TOW-2A  donors  against  palletized  M107  acceptors  are  shown  in  Figure  13. 
Neither  detonation  nor  burning  propagation  is  predicted.  The  other  probabilities  remain  high  and  then 
drop  rapidly  to  zero  with  increasing  range.  Probabilities  of  mechanical  damage  exceeding  1%  persist  to 
127  ft  The  1%  distance  for  hit  probability  is  662  ft. 

Similar  results  for  palletized  TOW-2A  donors  against  palletized  TOW-2A  acceptors  are  shown  in 
Figure  14.  No  detonation  propagation  is  predicted,  but  burning  propagation  occurs  at  the  1%  probability 
level  out  to  a  range  of  129  ft  Probabilities  of  mechanical  damage  of  1%  or  greater  persist  to  149  ft.  The 
1%  distance  for  hit  probability  is  667  ft  The  fact  that  the  bum  and  mechanical  damage  curves  nearly 
coincide  indicates  that  the  arbitrary  kinetic  energy  threshold  for  mechanical  damage  is  too  high. 

The  1%  propagation  probability  ranges  discussed  in  the  foregoing  paragraphs  are  summarized  in 
Table  4. 

The  lethality  of  the  three  donors  varies  substantially.  Significant  differences  between  the  single  and 
palletized  M107  donors  occur  simply  because  of  the  number  of  munitions  involved.  The  probabilities  for 
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Table  4.  1%  Propagation  Probability  Distances  (ft) 


Donor 

Acceptor 

Detonation 

Burning 

Mechanical 

Damage 

Hit 

Single  M107 

M107  Stack 

687 

mm 

TOW-2A  Stack 

iKm 

666 

HI 

M107  Stack 

M107  Stack 

642 

mm 

2,950 

TOW-2A  Stack 

643 

■sa 

TOW-2A  Stack 

M107  Stack 

NP 

NP 

mm 

662 

TOW-2A  Stack 

NP 

129 

Mm 

667 

NP  =  no  propagation 


these  donors  drop  slowly  and  the  maximum  ranges  are  considerably  greater  than  the  1%  ranges.  TOW-2  A 
donors  are  much  less  lethal  due  to  the  small  fragments  they  produce.  The  character  of  the  probability 
curves  for  these  donors  is  different.  Probabilities  remain  high  until  approaching  the  maximum  range 
where  they  drop  sharply. 

The  relative  vulnerability  of  the  two  acceptors  depends  on  the  type  of  response.  The  TOW-2A  is 
considerably  more  susceptible  to  burning  than  the  M107.  For  the  other  responses,  the  1%  propagation 
distances  are  generally  very  similar.  Although  the  (rather  arbitrary)  kinetic  energy  threshold  is  lower  for 
TOW-2A  acceptors,  their  probability  of  mechanical  damage  is  reduced  due  to  the  presence  of  the 
container. 

It  is  notable  that  this  analysis  (which  includes  a  number  of  worst-case  assumptions)  does  not  preclude 
propagation  of  detonation  from  M107  stacks  stored  at  the  required  253-ft  distance  from  either  acceptor. 
At  250  ft,  the  probability  of  detonation  propagation  to  an  M107  stack  is  between  68%  and  94% 
(depending  on  the  stack  configurations)  and  the  probability  of  detonation  propagation  to  a  TOW-2A  stack 
is  between  71%  and  91%. 

10.  SUMMARY  AND  CONCLUSIONS 

By  combining  several  existing  models,  we  have  developed  a  tool  (FRAGPROP)  for  estimating  the 
probabilities  associated  with  the  propagation  of  reaction  between  user-described  ammunition  stacks.  The 


27 


acceptor  responses  considered  include  detonation  of  energetic  materials,  burning  of  energetic  materials  (but 
not  of  combustible  packaging),  and  mechanical  damage.  The  models  include  the  FRAGHAZ  program  for 
the  Monte  Carlo  treatment  of  fragment  trajectories  and  the  accumulation  of  hit  probabilities,  the  Jacobs- 
Roslund  criterion  for  initiation  of  detonation,  and  the  ballistic  limit  condition  for  initiation  of  burning. 

A  number  of  difficulties  arose  in  applying  this  tool  to  munitions  of  interest.  Since  the  appropriate 
fragmentation  input  data  was  not  always  available,  notably  in  the  case  of  missiles,  we  developed  methods 
of  estimating  this  data.  It  was  necessary  to  represent  the  vulnerability  of  explosives  and  propellants  for 
which  Jacobs-Roslund  parameters  are  not  available  using  estimates  based  on  values  for  similar 
compositions.  It  was  also  necessary  to  assess  the  lethality  of  aluminum  fragments  with  models  intended 
for  use  with  steel  fragments. 

Although  only  two  weapons  have  been  specifically  analyzed,  the  responses  of  the  thick-walled  M107 
and  the  thin-walled  TOW-2A  are  expected  to  encompass  the  range  of  responses  of  a  wide  variety  of 
munitions.  The  analysis  was  used  to  predict  distances  below  which  the  propagation  probabilities  for  each 
response  exceed  1%  (see  Table  4).  Based  on  these  predictions,  we  conclude  that  the  artillety  ammunition 
donor  stacks  are  much  more  lethal  than  the  missile  donor  stacks  (which  cannot  produce  detonation)  and 
the  missile  acceptor  stacks  are  more  vulnerable  to  the  propagation  of  burning  (but  not  of  detonation)  than 
the  artillery  ammunition  acceptor  stacks.  However,  different  assumptions  regarding  the  orientation  of 
fragments  at  impact  time  might  enhance  sensitivity  to  burning. 
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APPENDIX  A: 

PROPAGATION  ACCIDENTS  IN  AMMUNITION  HOLDING  AREAS 
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We  considered  several  sources  of  information  which  are  pertinent  to  propagation  accidents  in 
ammunition  holding  areas.  These  include  reports  of  actual  events  found  in  the  files  of  the  Department 
of  Defense  Explosive  Safety  Board  as  well  as  results  of  hazard  classification,  packaging  and  MLLVAN 
tests. 


At  least  57  accidents  involving  open  storage  of  ammunition  have  been  identified.  We  have  reviewed 
reports  of  many  of  these.  There  are  a  number  of  difficulties  associated  with  the  data  given  in  these 
reports.  Much  of  it  is  old,  dating  to  the  World  War  n  era.  The  observations  were  made  at  different  times 
by  different  people  with  different  ideas  about  what  is  important.  The  amount  of  information  reported 
varies  widely,  and  details  of  storage  configurations  were  not  usually  included. 

We  drew  a  number  of  conclusions  from  our  review.  En  masse  detonation  of  large  quantities  of 
ammunition  is  rare,  usually  involving  bombs  or  large  projectiles.  More  commonly,  propagation  accidents 
follow  a  sequence  of  events  involving  fire  in  an  ammunition  stack  and  cookoff  of  munitions  in  that  stack. 
This  leads  to  propagation  of  fire  to  neighboring  stacks  via  hot  fragments  and  burning  debris.  The  process 
is  then  repeated  in  events  which  may  take  days  to  unfold.  Propulsive  reactions  in  which  items  such  as 
rockets  and  mortar  rounds  are  launched  toward  neighboring  or  distant  stacks  often  contribute  to 
propagation.  White  phosphorus  rounds  are  also  major  contributors  to  this  type  of  propagation.  Three 
illustrative  accident  reports  are  summarized  in  Table  A-l. 

The  most  relevant  hazard  classification  test  is  the  bonfire  test.  While  this  provides  the  maximum 
fragment  radius,  it  doesn’t  always  identify  the  number  of  fragments  which  travel  shorter  distances. 
Distances  over  which  firebrands  are  spread  are  usually  not  reported.  Another  relevant  series  of  tests  were 
conducted  in  MILVAN  containers.1  Although  the  quantities  of  explosive  in  each  MILVAN  were  modest, 
items  were  thrown  large  distances. 

The  importance  of  packaging  was  demonstrated  by  Teitell  and  Reeves2  who  conducted  tests  to 
determine  the  vulnerability  of  several  munitions  in  wood  packaging  to  several  threats.  The  general  results 


1  Lawrence,  W.  "Fragment  Hazards  From  Munitions  in  Containers."  BRL-TR-3203,  U.S.  Army  Ballistic  Research  Laboratory, 
Aberdeen  Proving  Ground,  MD,  1991. 

2 

Teitell,  L.,  and  H.  J.  Reeves.  ’’Fire  Retardant  Packaging  for  Artillery  Ammuniton."  BRL-MR-2490,  U.S.  Army  Ballistic 
Research  Laboratoiy,  Aberdeen  Proving  Ground,  MD,  1975. 
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Table  A-l.  Summary  of  Representative  Propagation  Accident  Reports 


DDESB  FILE 

466 

480 

1271 

LOCATION 

PUEBLO,  CO 

IESHIMA,  JAPAN 

DA  NANG,  VIETNAM 

AMMUNITION 

3  "PILES"  OF 
TNT-LOADED  3" 
M42  AMMUNITION 
SEPARATED  BY 
120'  AND  285' 

2162  TONS  OF 
BOMBS,  ROCKETS, 
SMALL  ARMS,  20mm 
IN  OPEN  STORAGE 
WITH  ABOUT  200’ 
BETWEEN  STACKS 

LARGE  MARINE 
CORPS  AMMUNITION 
SUPPLY  POINT 

CAUSE 

LIGHTNING 

FIRE,  POSSIBLY 
CAUSED  BY 
ENEMY  ACTION 

DURATION 

HOURS 

2  DAYS 

>  24  HOURS 

RESULT 

•MULTIPLE  EVENTS 

•  FIRE  SPREAD  BY 
HOT  FRAGMENTS 

AND  BURNING 
DEBRIS  *VIRTUALY 
ALL  AMMO  LOST 

•  FRAGS  TO  3700' 

•  THROWN  ROUND 
DETONATED  ON 
IMPACT  AT  200‘ 

•  BURNING  WOOD 
AND  FIBERBOARD 

TRAVELED  400'. 

ESSENTIALLY 
EVERYTHING  LOST 

•  MULTIPLE 
EXPLOSIONS 

•  ROCKETS 
COMMUNICATED 
REACTION  TO  AN 
AIR  FORCE  ASP 

•  LARGE  LOSS  OF 
AMMO  •  "NO  CELL 
APPEARED  TO 
COMMUNICATE 
DIRECTLY  TO  ITS 
NEIGHBOR." 

of  these  tests  were  that  the  threats  did  not  cause  detonations  but  started  fires  in  propellants  and  wood 


packaging.  Propellant  fires  rapidly  cooked  off  warheads  and  wood  fires  cooked  off  warheads  causing  both 
detonations  and  less  violent  explosions.  Detonations  sometimes  scattered  stacks  and  stopped  reaction. 
Fire-resistant  packaging  prevented  the  spread  of  reaction  in  some  cases. 


We  conclude  that  the  dominant  mode  of  propagation  in  ammunition  accidents  is  fire  leading  to  cookoff 
leading  to  more  fire  and  more  cookoff.  Combustible  packaging  and  propulsive  reactions  are  major 
contributors  to  propagation.  Insensitive  munition  technology  will  reduce  the  probability  of  propagation, 
especially  where  stacks  can  be  separated  by  more  that  50  ft.  However,  quantitative  analysis  of  propagation 
probability  for  this  mechanism  is  probably  not  possible  at  this  time. 
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APPENDIX  B: 
FRAGPROP  LISTING 
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ooo  oooooo  oo  oo  ooo  oooo  ooooooooo  oo  oooono 


*********************************************************************** 
program  frgpr 


FRAGPROP 


common  blocks 

common/cnstm/  qpi,  hpi,pi, tpi,  raddeg,  degrade  smlrad,  smlnum, bignum 
common/ repls/  nreps,fnreps 
common/frags/  nf rags, tf rags 

common/ range/  rngmax,  rngseg,  rngsgh,  rngscl,  rngsix,  rngmin 
common/ rkvar/  rkv(7) 

type  declarations 
character* 8  unielaf unirem 
double  precision  rndinp, rndnxt, rnds a v 

report  the  program  title 
write (*, '  (/, a) ' )  '  FRAGPROP' 

set  up  the  run 

call  setup (rndinp, rndsav, rndnxt,  ihbgn,tsbgn) 

begin  the  replication  loop 
do  300  irep=l,nreps 

if rcm=  (irep-1) *nfrags 

initialize  the  replication  parameters 

call  rpini (rndnxt,  slcnst, asldegf alttud, wspeed, wdrdeg, wdrrad, 
+  jxrng, jrkmax) 

begin  the  fragment  loop 
do  200  ifrag=l, nf rags 

estimate  the  remaining  time 

frcmp=max (float (if rcm+if rag-1) /tf  rags,  smlnum) 

call  ertim (ihbgn,  tsbgn,  f rcmp,  timela, uniela, timrem,  unirem) 

report  the  current  replication  and 
fragment  and  the  remaining  time 
write (*, ' (a,i3,a,i3,a,i3,a,i3,a,f5.2,a)') 

+  f +Replication:  9  f irepf r /' , nrepsf 9  Fragment:  f,ifragr'/ff 
+  nfragsff  Remaining  Time:  9  /timrem, unirem 

initialize  the  fragment  parameters 

call  frini (irepf ifragf  rndnxt ,  irico, lrangef htf rft ,mf rhtlf 
+  vf ragf  seldeg, selrad, aeldeg, aelradf  cdlf  cd2, cd3f 

+  deldst) 

begin  the  trajectory  computation  loop 
itraj=0 

100  itra j=itra j+1 

determine  the  time  step 

call  tstep (deldst f iricof vfragf aelrad,  deltim) 

increment  the  time 
rkv ( 1 ) =rkv ( 1 ) +deltim 

perform  the  Runge-Kutta  integration 

call  rkint (wspeed, wdrdeg, wdrrad, jxrng, jrkmax, cdl, cd2, cd3, 
+  deltim, alttud,  if rag) 

determine  the  fragment  variable  values 
at  the  end  of  the  integration  step 

call  f rvar (vf rag, selrad, seldeg, aelrad, aeldeg,  range, irange) 


F000010 

F000020 

F000030 

F000040 

F000050 

F000060 

F000070 

F000080 

F000090 

F000100 

F000110 

F000120 

F000130 

F000140 

F000150 

F000160 

F000170 

F000180 

F000190 

F000200 

F000210 

F000220 

F000230 

F000240 

F000250 

F000260 

F000270 

F000280 

F000290 

F000300 

F000310 

F000320 

F000330 

F000340 

F000350 

F000360 

F000370 

F000380 

F000390 

F000400 

F000410 

F000420 

F000430 

F000440 

F000450 

F000460 

F000470 

F000480 

F000490 

F000500 

F000510 

F000520 

F000530 

F000540 

F000550 

F000560 

F000570 

F000580 

F000590 

F000600 

F000610 

F000620 

F000630 

F000640 

F000650 

F000660 

F000670 

F000680 
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oo  o  o  oo  oo  oo  oo  oo  oo  oo  o  oo  oo  oo  oo  oo  oo  oo  oo 


perform  hazard  computations  as  required 

call  hazch  (if rag,  irange,mfrhtl,  mfrht2,  lrange,  aelrad,  vfrag) 

check  the  range  for  trajectory  termination 
if (range. gt.rngmax)  go  to  200 

check  the  fragment  velocity  for  trajectory  termination 
if (vfrag. It. 20.0)  go  to  200 

check  for  ricochet 
if (rkv(5) .le.0.0)  then 

check  the  ricochet  count 
if (irico.gt. 6)  go  to  200 

check  the  elevation  angle  for  trajectory  termination 
if (aeldeg.ge . asldeg)  go  to  200 

determine  the  ricochet  velocity  and  angle 

call  ricoc (slcnst, vfrag, aeldeg, selrad, vrico,  ricrad) 

check  the  ricochet  velocity  for  trajectory  termination 
if (vrico. It. 20.0)  go  to  200 

reinitialize  the  velocity  integration  variables 
call  reini (irico, vrico, ricrad, wspeed, wdrrad,  vfrag, 

+  selrad, seldeg, aelrad, aeldeg, deldst,mf rhtl) 

endif 

end  of  the  trajectory  computation  loop 
go  to  100 

end  of  the  fragment  loop 
200  continue 

accumulate  hit  probabilities 
call  hpacc 

end  of  the  replication  loop 
300  continue 

determine  the  elapsed  time 

call  ertim (ihbgn, tsbgn, 1.0, timela, uniela, timrem, unirem) 

write  the  output  files 
call  outpt 

report  completion 
write (*, ' (a,f5.2,a/) ' ) 

+  ' +Computation  Complete  Total  Elapsed  Time :  ' , 

+  timela, uniela 

stop  '  ' 

end 
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Q** ********************************************************** *********** 

c 

subroutine  setup (rndinp, rndsav, rndnxt ,  ihbgn, tsbgn) 
c 

c  sets  up  the  run 

c 

c  common  blocks 

common/error/  meror, irerr, iferr, merout 
common/rangs/  nrsegs, nrsegp 
common/denst/  fd(97) ,fdd(97) f fdb(97) , fdm(97) 
common/phitt/  ph(97) ,phd(97) ,phb(97) ,phm(97) 

type  declarations 
character  dfile*25f afile*25 
double  precision  rndinp,  rndnxt,  rndsav 

assign  constants 
call  const 

obtain  input 

call  input (rndinp, rndsav, rndnxt, dfile, minvrt, afile) 

compute  parameters  requiring  both  donor  and  acceptor  input 
call  doacc 

output  the  run  conditions 
call  out in (rndsav, df ile, minvrt , afile) 

compute  the  maximum  vulnerable  areas  of  the  acceptor 
call  vulna 

initialize  the  hit  probabilities 
call  arini (nrsegs, ph,phd,phb,phm, 0.0) 

initialize  the  fragment  densities 
call  arini (nrsegs, fd, fdd, fdb, fdm, 0.0) 

obtain  the  starting  time 
call  gettim ( ihbgn, imbgn, isbgn, ilbgn) 
tsbgn=3600 . 0*f loat (ihbgn) +60 . 0*f loat (imbgn) 

+  +f loat  (isbgn) +f loat (ilbgn) /100.0 
c 

c  report  return 

if (meror .gt . 0)  write (2, '  (lx, a)')  'returning  from  setup' 


end 

c 

c** ************** ****************************************** ************* 
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c 

subroutine  const 
c 
c 
c 
c 


c 

c 


c 

c 


c 

c 


c 

c 


assigns  constants 
common  blocks 


common/cnstm/  qpi , hpi , pi , tpi,  raddeg, degrad, smlrad,  smlnum, bignum 
common/ cnstp/  gg, tglbgr , hrair 0 , racnst , vsound, hccnst , hcaml , hcam2 , 
+  hcam3 

common/ cnstu/  fgrtlb, fgmtlb, fgrtgm, ffttm, ffttmm, fmtft, f i2tc2 

mathematical  constants 
qpi=atan<1.0) 
hpi=2 . 0  *qpi 
pi=4 . 0*qpi 
tpi=2 .  0*pi 
raddeg=360 . 0/tpi 
degrad=tpi/360 . 0 
smlrad=l . Oe-3  *degrad 
bignum=l . 0e30 
smlnum=l . 0e-30 

physical  constants 
gg=32.174 
tglbgr=l . 4e4*gg 
hrair0=3 . 8239e-2 
racnst=3 . 39e4 
vsound=l . 1164e3 
hccnst=2 . 86e5 
hcaml=0.75 
hcam2=l . 5 
hcam3=2 . 5 

unit  conversions 
fgrtlb=l./7.0e3 
f gmtlb=2 . 2046226e-3 
fgrtgm=fgrtlb/fgmtlb 
ffttm-0.3048 
f fttmm=l . e-3*f fttm 
fmtft=l . 0/f fttm 
fi2tc2=6 .4516 

return 

end 
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c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  input (rndinpf rndsav, rndnxt , dfile, minvrtf af ile) 

obtains  input 

type  declarations 
character  dfile*25f af ile*25 
double  precision  rndinpr rndnxt, rndsav 

input  randomization  seed 
call  rndin (rndinp, rndsav, rndnxt ) 

input  parameters 
call  parin 

input  donor  data 
call  donin (dfile, minvrt) 

input  acceptor  data 
call  accin(afile) 

input  range  data 
call  rngin 


return 


c 


end 
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***********************************************************************  Rqoooio 


R000020 

subroutine  rndin (rndinp, rndsav, rndnxt)  R000030 

R000040 

obtains  randomization  seed  input  R000050 

R000060 

type  declarations  R000070 

double  precision  rndinp, rndsav, rndnxt  R000080 

R000090 

100  write (*,' (a\)')  '  Enter  Monte  Carlo  Seed  (1  to  2147483646):  '  R000100 

read{*, *,err=100)  rndinp  R000110 

if (rndinp. It. 1.0. or. rndinp. gt. 2147483646.0)  go  to  100  R000120 

R000130 

rndsav=rndinp  R000140 

rnddum=rndom (rndinp)  R000150 

rndnxt=0 . OdO  R000160 

R000170 

return  R000180 

R000190 

end  R000200 

R000210 


**************************************************** *****************  ro 00220 


c 

c 

c 

c 


c 


c 

c 

c 


*********************************************************************  P000010 


P000020 

subroutine  parin  P000030 

P000040 

obtains  parameter  input  P000050 

P000060 

common  blocks  P000070 

common/error/  meror, irerr, iferr,merout  P000080 

common/repls/  nreps,fnreps  P000090 

common/prmtr/  slcmin, slcmax, altmin, altmax,  P000100 

+  wspmin, wspmax, wdrmin, wdrmax  P000110 

P000120 

open(4,file=' fp.prm' , status=' old' )  P000130 

read(4,*)  meror, irerr, iferr  P000140 

if (meror. gt.0)  open (2, file=' fp. err' , status=' unknown' )  P000150 

read (4, *)  nreps  P000160 

nreps=min (nreps, 999)  P000170 

fnreps=float (nreps)  P000180 

read(4,*)  slcmin, slcmax  P000190 

read(4,*)  altmin, altmax  P000200 

read (4,*)  wspmin, wspmax  P000210 

read(4,*)  wdrmin, wdrmax  P000220 

close (4)  P000230 

P000240 

return  P000250 

P000260 

end  P000270 


P000280 


c***********************************************************************  poo 0290 
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c** ************** **************************** ****************** *********1)000010 


c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 


c 

c 


subroutine  donin (df ilef minvrt) 

obtains  donor  input 

common  blocks 

common /cnstm/  qpi, hpif  pi,  tpi,  raddeg,  degrad,  smlrad,  smlnum, bignum 
common/cnstp/  gg, tglbgr, hrairO, racnst,  vsound,hccnst,hcaml,hcam2, 
+  hcam3 

common/ cnstu/  fgrtlb,  fgmtlb,  fgrtgm,  ffttm,  ffttmm,  fmtft, f i2tc2 

common/error/  meror, irerr, iferr,merout 

common/f  rmin/  fmgrmn,  fmgmmn,  ezndeg 

common/zones/  elzdegf  azsdeg,  azsrad 

common/donor/  nunits,  htstft f  htbaft 

common/frags/  nf rags, tf rags 

common/f ragd/  famil{900) , vfifs  (900) , farma(900) , 

+  fatn (900) , ezldeg(900) 

common/fragc/  facm2  (900) ,  fdmm(900) ,  fmgm(900) 
common/f ragx/  vcrdx(900) , vcrbx(900) 
common/ f ragk/  frkec (900) 
common/names/  dname,aname 

local  array 
dimension  fmult(36) 

type  declarations 
character  dname*40, aname*40 
character  dfile*25, answr*4 

report  call 

if (meror. gt.0)  write (2, ' (4x, a) ' )  'donin  called' 

user  specification  input 
open (4, f ile=' fp.dnr' , status=' old' ) 
read (4, ' (40x,a) ' )  dname 
read (4, '  (4 0xr a) ' )  dfile 
read(4f ' (40x, a) ' )  answr 
minvrt=0 

if (answr . eq. ' Y' )  minvrt=l 
if (answr .eq. 'y' )  minvrt=l 
if (answr .eq. 'YES' )  minvrt=l 
if (answr .eq. 'Yes' )  minvrt^l 
if (answr .eq. 'yes' )  minvrt=l 
read(4,' (bnf 40xf i5) ' )  nunits 
read(4,' (bn, 40x,fl0.2) ' )  htstft 
read(4,' (bnf 40x, f 10 .2) ' )  htbaft 
read (4, ' (bn, 40x, fl0.2) ' )  fmgrmn 
fmgmmn=fgrtgm*fmgrmn 
read(4, ' (bn, 40x, flO .2) ' )  ezndeg 
close  (4) 

donor  lethality  input 
open (4, file=dfile, status='old' ) 
read ( 4 , * )  elzdeg, azsdeg 
azsrad=degrad*azsdeg 
read (4, *)  nfmlts,mdosh 
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if (mdosh.gt • 0)  then 
read (4,*)  dshin 
dshcm-2 . 54*dshin 

read {4, *)  ccdsv, athdsv,bthdsv,  cthdsv 
tcdsv=10 . 0**ccdsv 

read ( 4 , * )  ccdsm, athdsm, bthdsm, cthdsm 
tcdsm=10 . 0**ccdsm 
endif 

read(4, *)  (fmult (ifmlt) , ifmlt=l,nfmlts) 
fnuaz=f loat (nunits) *azsdeg 
if rag=l 

100  read(4, *, end=200)  poldeg,  fmgr,tf  amil,  tvf  if  srtf  arma 
if (minvrt .eq. 0)  then 
tezdeg=90 . 0-poldeg 
else 

tezdeg=poldeg-90 . 0 
endif 

if (tezdeg.lt .ezndeg)  go  to  100 
tfmgm=0 . 0647 9891 *fmgr 
tfacm2=f i2tc2*tfamil*fgrtlb*fmgr 
if (mdosh . gt . 0 )  then 

call  vmthr  (dshcm,  tf acm2 , tvf if s,  tfmgm,  teds v,  athdsv, bthdsv,  cthdsv, 
+  tcdsmf  athdsmf  bthdsmr  cthdsmf  tvf  if  sf  tfmgm) 

if (tvf if s . le . 0 . 0)  go  to  100 
endif 

if  (tfmgm.lt  .fmgmmn)  go  to  100 
if  (ifrag.gt . 900)  then 
write (*, ' (a, i3, a) 9  ) 

+  9  Data  file  contains  more  than  ',900,'  fragments!' 

stop  'Computation  Terminated' 
endif 

ezldeg (if rag) =tezdeg 
fmgm(if  rag)  =tfmgm 
f amil (if rag) =tf amil 
vf if s (if rag) =tvf if s 
f  arma  (if  rag)  =tf  arma 

fatn(ifrag)=fnuaz*fmult (int (poldeg/elzdeg) ) 
f acm2 (if rag) =tf acm2 

fdmm(if rag) =2 . 0*sqrt  (100 . 0*tfacm2/pi) 
f rkec (if rag) =fmgr/tglbgr 
ifrag=ifrag+l 
go  to  100 

200  nfrags=if rag-1 
close (4) 

report  return 

if (meror .gt . 0)  write (2f' (4x/a)')  'returning  from  donin' 
return 


end 

c 

q** ********************************************************** *********** 


c 

c 
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c****** ******************************************************** *********  AO 000 10 
c  A000020 


subroutine  accin(afile) 

A000030 

c 

A000040 

c 

obtains  acceptor  input 

A000050 

c 

A000060 

c 

common  blocks 

A000070 

common/error/  meror, irerr, if err f merout 

A000080 

common/names/  dname,  aname 

A000090 

common/acptr/  htacft , wdacf t , dpacf t , htacfm 

A000100 

common/ acf lg/  nweps , macwh,macrmf macsh , mstor 

A000110 

common/geowp/  wplin,wpdin 

A000120 

common/geopa/  pasin,pahin, pasft 

A000130 

common / geo wh/  whlin,  whdin,  whhin,  whhcm,  whhmm,  exlin 

A000140 

common/ jrdtx/  abjrx,  cjrx 

A000150 

common/ thorw/  athw,bthw,  cthwf  dthw,  vrwms,  tew 

A000160 

common/georm/  rmlin,  rmdin,  rmhin,  rmhem,  rmhmm,  prlin 

A000170 

common/ jrdtp/  abjrp,cjrp 

A0  00 180 

common/ thorr/  athr , bthr , cthr,  dthr, vrrms , ter 

A000190 

common/geoas/  ashem 

A000200 

common/ thras/  athasv, bthasv, cthas v,  tcasv. 

A000210 

+  athasm,  bthasmf  cthasm,  tcasm 

A000220 

common /mechd/  wpke 

A000230 

c 

A000240 

c 

type  declarations 

A000250 

character  dname*40f aname*40 

A000260 

character  afile*25 

A000270 

c 

A000280 

c 

report  call 

A000290 

if (meror .gt . 0)  write (2f r  (4xf a) ' )  'accin  called' 

A000300 

c 

A000310 

c 

user  specification  input 

A000320 

open (4f  f ile=' fp.acc' , status^' old' ) 

A000330 

read (4, ' <40x, a) ' )  aname 

A000340 

read (4, ' (40xf a) ' )  afile 

A000350 

read(4f ' (bn, 40x, f 10 .2) ' )  htacft 

A000360 

htacfm=htacf t~0 . 2 

A0  00 370 

read(4f '  (bnf  40xf  flO .2) ' )  wdaeft 

A000380 

read(4f ' (bnf 40xf f 10 .2) ' )  dpaeft 

A000390 

close (4) 

A000400 

c 

A000410 

c 

acceptor  vulnerability  input 

A000420 

open (4f  f ile=af iler  status=' old' ) 

A000430 

read  ( 4 ,  * )  nweps , macwh, maerm, macsh, mstor 

A000440 

read(4f  *)  wplin, wpdin 

A000450 

read ( 4  f  * )  pasin, pahin 

A000460 

pasft=pasin/12 . 

A0  00 470 

if (macwh . ne • 0 )  then 

A000480 

read (4, *)  whlin, whdin, whhin, exlin 

A000490 

whhcm=2 . 54*whhin 

A000500 

whhmm=10 . 0*whhcm 

A000510 

read (4, *)  a jrx,b jrxpl, c jrx 

A000520 

ab  j  r x=a  j  r x  *b  j  r xp 1 

A000530 

read  ( 4 , * )  ccw, athw, bthw, cthw, dthw, vrwms 

A000540 

tew- 10 . 0**ccw 

A000550 

endif 

A000560 

if (macrm.ne.0)  then 

A0 00 570 

read (4,  *)  rmlin,  rmdin,  rmhin, prlin 

A000580 

rmhcm=2 . 54*rmhin 

A000590 

rmhmm=10 . 0*rmhcm 

A000600 

read(4, *)  a jrp,bjrppl, c jrp 

A000610 

abjrp=a jrp*b jrppl 

A000620 
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read (4, *)  ccr f athr, bthr, cthr , dthr, vrrms 
tcr=10 . 0**ccr 
endif 

if (macsh . ne . 0 )  then 
read (4,*)  ashin 
ashcm=2 . 54*ashin 

read (4, *)  ccasv, athasv,bthasv, cthasv 
tcasv=10 . 0**ccasv 

read ( 4 , * )  ccasm,  athasm, bthasm,  cthasm 
tcasm=10 . 0**ccasm 
endif 

read ( 4 , * )  wpke 
close  (4) 
c 

c  report  return 

if (meror .gt . 0)  write (2, ' (4x, a) ' )  'returning  from  accin' 


end 

c 

c** ******************************************************** ************* 
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c** ********************************************************** *********** 
c 

subroutine  rngin 
c 

c  obtains  range  input 

c 

c  common  blocks 

common/error/  meror, irerr, if err ,  merout 
common/rangs/  nrsegs,nrsegp 

common/ range/  rngmax, rngseg, rngsgh, rngscl, rngsix, rngmin 
common/arang/  rangem(97) ,ahfrn<97) ,ahtop(97) 
c 

c  report  call 

if (meror. gt.0)  write (2, '  (2x, a) ' )  f rngin  called' 

user  specification  input 
open ( 4 , f ile=f  fp . rng' , status=f  old' ) 
read (4, ' (bn, 4 Ox, f 10 .2) ' )  rngmax 
read (4, ' (bn, 40x, i5) ' )  nrsegs 
nrsegs=min (nrsegs, 96) 
nrsegp-nrsegs+1 
rngseg=rngmax/f loat (nrsegs) 
rngsgh=0 . 5*rngseg 
rngscl=l . 0e-2*rngseg 
rngsix=0 . 6*rngseg 
do  100  irange=l,nrsegp 

rangem (i range) ^rngseg* (float (irange) -0 . 5) 

100  continue 
close (4) 

report  return 

if (meror. gt.0)  write (2, '  (4x, a) ' )  'returning  from  rngin' 
return 
end 
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subroutine  doacc 

computes  parameters  requiring  both  donor  and  acceptor  input 
common  blocks 

common/error/  meror, irerr, if err f merout 
common/repls/  nreps,fnreps 
common/ zones/  elzdeg,  azsdeg, azsrad 
common/frags/  nfrags, tf rags 

common/fragd/  famil (900) , vf ifs (900) , farma (900) , 

+  fatn(900) f ezldeg(900) 

common/fragc/  facm2 (900) , fdmm (900) , fmgm{900) 
common/ fragx/  vcrdx (900) , vcrbx (900) 
common/fragp/  vcrdp (900) , vcrbp (900) 
common/acf I g/  nweps , macwh, macrm, macsh, mstor 
common/acptr/  htacft, wdacft ,dpacft. htacfm 
common/geowh/  whlin,  whdin,  whhm,  whhcm,  whhmm,  exlin 
common/ jrdtx/  abjrx, cjrx 

common/ thorw/  atnw^btnw^  cthw, dthw, vrwms, tew 
common/ georm/  rmlin, rmdin, rmnin, rmhem, rmhmm, prlin 
common/ thorr/  athr, bthr, cthr, dthr , vrrms, ter 
common/ jrdtp/  abjrp,cjrp 
common/rangs/  nrsegs, nrsegp 

common/ range/  rngmax, rngseg, rngsgh, rngscl , rngsix, rngmin 
common/arang/  rangem(97) , ahfrn{97) , ahtop(97) 

report  call 

if (meror.gt.0)  write (2, '  (2x, a) ' )  'doacc  called' 

compute  the  total  number  of  fragment  trajectories 
tfrags=float (nreps*nfrags) 

compute  the  minimum  range  and  range  area  arrays 
tahazs=tan (0 . 5*azsrad) 
if (tahazs .gt . 0 . 0)  then 

rngmin=wdacft/ (2.0*tahazs) 
else 

rngmin=0 . 0 
endif 

do  100  irange=l, nrsegp 

rangeo=rangem(irange) + rngsgh 

rangei=rangem (irange) -rngsgh 

ahfrn (irange) =azsrad*htacft*rangem {irange) 

ahtop (irange) =0 .5*azsrad* (rangeo*rangeo-rangei*rangei) 

100  continue 

?  redetermine  parameters  characterizing 
ragment  lethality  against  unshrouded  acceptor 
if (macsh. eq. 0)  then 

if {macwh. ne . 0)  then 

a  warhead  component  is  present 
do  200  ifrag=l, nfrags 

vcrdx (if rag) =vcr jr (ab jrx, cjrx,  whhmm, fdmm (if rag) ) 
vcrbx (if rag) = 

+  vcrth(tcw, athw, bthw, whhcm, facm2 (if rag) , fmgm(ifrag) ) 

200  continue 

endif 

if (macrm. ne .0)  then 

a  rocket  motor  component  is  present 
do  300  ifrag=l, nfrags 

vcrdp (ifrag)=vcr jr (abjrp,  cjrp, whhmm, fdmm(ifrag) ) 
vcrbp (if rag) = 

+  verth (ter, athr,bthr, rmhem, facm2 (if rag) , fmgm(ifrag) ) 

300  continue 

endif 

endif 


report  return 

if (meror.gt.0)  write (2, ' (2x, a) ' )  'returning  from  doacc' 

return 

end 
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***********************************************************************  V000010 


subroutine  vulna 

computes  the  maximum  vulnerable  areas  of  the  acceptor 
common  blocks 

common/ error/  meror , irerr , if err,  merout 
common/ acflg/  nweps,macwh,macrm,macsh,mstor 
common/ vulns/  amf rn, amtop, axf rn, axtop, apf rn, aptop 

report  call 

if (meror. gt.O)  write (2, '  (2x, a) ' )  'vulna  called' 

initialize  the  vulnerable  areas 
amf rn=0 . 0 
amtop=0 . 0 
axf rn=0 . 0 
axtop=0 . 0 
apf rn=0 . 0 
aptop=0 . 0 

if (mstor.eq. 1)  then 

the  storage  arrangement  is  vertical 
call  vulnv(ucfrn,uctop) 

if  a  warhead  component  is  present 
if (macwh.ne.O)  call  vulvw(ucfrn,uctop) 

if  a  rocket  motor  component  is  present 
if (macrm.ne.O)  call  vulvm(ucfrn,uctop) 

else 

the  storage  arrangement  is  horizontal 
call  vulnh (ucf rn,uctop) 


amin2=0 . 0 

if  a  warhead  component  is  present 
if (macwh.ne.O)  call  vulhw(amin2,ucfrn,uctop) 

if  a  rocket  motor  component  is  present 
if (macrm.ne . 0)  call  vulhm(amin2, ucfrn,uctop) 

amf  rn=ucfrn*amin2 
amt  op=uct  op*amin2 

endif 

report  return 

if (meror. gt.O)  write (2, ' (2x,a) ' )  ' returning  from  vulna' 

return 

end 
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*********************************************************************** 


subroutine  vulnv (ucfrn, uctop) 

determines  the  number  of  front  and  top 
units  for  a  vertical  storage  arrangement 

common  blocks 

coramon/error/  meror, irerr, if err f merout 
common/acptr/  htacf t , wdacf t , dpacf t , htacfm 
coinmon/geowp/  wplin,  wpdin 
common/geopa/  pasin,pahin,pasft 

report  call 

if (meror. gt .0)  write (2, ' (4x, a) '  )  'vulnv  called' 

compute  the  front  and  top  unit  areas  in  sq  in 
wdpps=wpdin+pasin 
auf rn= (wplin+pahin) *wdpps 
aut op=wdpps  *wdpps 

compute  the  numbers  of  front  and 
top  units  (with  conversion  to  sq  ft) 
wdacp=wdacft+pasft 
ucfrn=wdacp*htacft/auf rn 
uctop=wdacp* (dpacft+pasft) /autop 

report  return 

if (meror. gt.O)  write (2, ' (4x, a) ' )  'returning  from  vulnv' 

return 

end 

****************************************************************** 
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subroutine  vulvw (ucfrn, uctop) 

determines  the  front  and  top  vulnerable 
warhead  areas  for  vertical  storage 

common  blocks 

common/error/  meror, irerr, if err, merout 

common/cnstm/  qpi, hpi, pi, tpi,  raddeg, degrad,  smlrad,  smlnum, bignum 
common/ geowh/  whlin, whdin, whhin, whhcm, whhmm, exlin 
common/vulns/  amf rn, amtop, axf rn, axtop, apf rn, aptop 

report  call 

if (meror. gt . 0)  write (2, ' (4x, a) ' )  'vulvw  called' 

compute  the  front  and  top  maximum 
vulnerable  warhead  areas  in  sq  ft 
amfrn=ucfrn* whlin* whdin 
amtop=uctop*qpi*whdin*whdin 

compute  the  front  and  top  maximum 
vulnerable  explosive  areas  in  sq  ft 
exdin=whdin-whhin 
axfrn-ucfrn*exlin*exdin 
axtop=uctop*qpi*exdin*exdin 

report  return 

if (meror. gt.O)  write (2,' (4x, a)')  'returning  from  vulvw' 

return 

end 

t* ********************************************************** *********** 
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*********************************************************************** 

subroutine  vulvm  (ucfrn,uctop) 

determines  the  front  and  top  vulnerable 
rocket  motor  areas  for  vertical  storage 

(there  are  no  top  vulnerable  areas  if  a  warhead  is  also  present) 
common  blocks 

common/cnstm/  qpi,hpi,pi, tpi,  raddeg, degrad,  smlrad,  smlnum, bignum 
common/error/  meror, irerr, iferr,merout 
common/georm/  rmlin,  rmdin,  rmhin,  rmhcm,  rmhmm,prlin 
common/ vulns/  amf rn, amtop, axf rn, axtop, apf rn, aptop 

report  call 

if (meror. gt.O)  write (2, '  (4x, a) ' )  'vulvm  called' 

compute  the  front  and  top  maximum 
vulnerable  motor  areas  in  sq  ft 
amf  rn=amf  rn+ucf  rn  *  rml in *  rmdin 
if (amtop.eq. 0 . 0)  amtop=uctop*qpi*rmdin*rmdin 

compute  the  front  and  top  maximum 
vulnerable  propellant  areas  in  sq  ft 
prdin=rmdin-rmhin 
apf  rn=uc  frn*prlin  *prdin 

if ( axtop. eq. 0 . 0)  aptop=uctop*qpi*prdin*prdin 


report  return 

if (meror .gt . 0)  write (2, ' (4x, a) ' )  'returning  from  vulvm' 

return 

end 
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***********************************************************************  V0003 60 


**************************************** ******************** *********** 


subroutine  vulnh (ucfrn,uctop) 

determines  the  number  of  front  and  top  units  for  a 
horizontal  storage  arrangement  (nweps  weapons  high) 

common  blocks 

common/error/  meror, irerr, iferr,merout 
common/acptr/  htacft, wdacft,dpacft,htacfm 
common/acf lg/  nweps, macwh,macrm,macsh,mstor 
common/geowp/  wplm,wpdin 
common/geopa/  pasin,pahin,pasft 

report  call 

if (meror .gt . 0)  write (2, ' (4x, a) ' )  'vulnh  called' 

compute  the  top  and  front  unit  areas  in  sq  in 
autop=wplin* (wpdin+pasin) 
aufrn=autop+wplin*pahin/float (nweps) 

compute  the  numbers  of  front  and 
top  units  (with  conversion  to  sq  ft) 
ucf rn=wdacft* (htacft+pasft) /aufrn 
uctop=wdacft* (dpacft+pasft) /autop 

report  return 

if (meror .gt . 0)  write (2, ' (4x, a) ' )  ' returning  from  vulnh' 


end 
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************************** ************************************** ******* 

subroutine  vulhw (amin2f  ucf rnf  uctop) 

determines  the  front  and  top  vulnerable 
warhead  areas  for  horizontal  storage 

common  blocks 

common/error/  meror, irerrf if err, merout 

common/ geowh/  whlin,  whdin, whhin, whhcm, whhmm,  exlin 

common/ vulns/  amf  rn,  amt  op,  axf  rn,  axtop,  apf  rn,  apt  op 

report  call 

if (meror.gt .0)  write (2, ' (4x, a) ' )  ' vulhw  called' 

compute  the  front  and  top  maximum 
vulnerable  warhead  areas  in  sq  ft 
amin2=whlin*whdin 

compute  the  front  and  top  maximum 
vulnerable  explosive  areas  in  sq  ft 
axin2= (whdin-whhin) *exlin 
axf  rn=uc  f  r n  *  axin2 
axtop=uctop*axin2 

report  return 

if (meror.gt . 0)  write (2, ' (4x, a) '  )  'returning  from  vulhw' 

return 

end 
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***********************************************************************  V000010 


subroutine  vulhm (amin2 , ucf rn, uctop) 

determines  the  front  and  top  vulnerable 
rocket  motor  areas  for  horizontal  storage 

common  blocks 

common/error/  meror, irerr, if err, merout 
common/georm/  rmlin,  rmdin,  rmhin,  rmhcm,  rmhmm,prlin 
common/ vulns/  amf rn, amtop, axf rn, axtop, apf rn, aptop 

report  call 

if (meror.gt. 0)  write(2,' (4x,a)')  'vulhm  called' 

compute  the  front  and  top  maximum 
vulnerable  motor  areas  in  sq  ft 
amin2 =amin2 +rml  in  *  rmdin 

compute  the  front  and  top  maximum 
vulnerable  propellant  areas  in  sq  ft 
apin2= (rmdin-rmhin) *prlin 
apf  r n=ap  f  r  n+uc  frn*apin2 
apt  op=apt  op+uc t op  *  apin2 

report  return 

if (meror.gt .0)  write (2,' (4x,a)')  'returning  from  vulhm' 

return 

end 


*********************************************************************** 
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************************************************** ********************* 


subroutine  out in (rndsav, dfile,  minvrt , afile) 
outputs  problem  input 
common  blocks 

common/repls/  nrepsf fnreps 
common/f rmin/  fmgrmn,  fmgmmn,  ezndeg 
common/zones/  elzdeg, azsdegf azsraa 
common/prmtr/  slcmin,  slcmax,  altmin,  altmax, 

+  wspmin,  wspmax,  wdrmin,  wdrmax 

common/donor/  nunits,htstft,htbaft 
common/frags/  nf rags, tfrags 
common/rangs/  nrsegs, nrsegp 

common/range/  rngmax, rngseg, rngsgh,  rngscl,  rngsixf  rngmin 

common /names/  dname,aname 

common/ acpt r /  htacf t , wdacf t , dpacf t , ht acf m 

type  declarations 
character  dname*40, aname*40 
character  dfile*25,afile*25 
double  precision  rndsav 

open (3f  file-' fp.out' , status=' unknown' ) 

write (3,' (//33x,a)')  ' ***FRAGPROP***' 
write ( 3 , ' (  /1 4x,  a ) ' ) 

+  'REACTION  PROPAGATION  PROGRAM  FOR  FRAGMENTING  MUNITIONS' 

write (3, ' (//33x, a/) ' )  ' **PARAMETERS**' 

call  inout ( '  NUMBER  OF  REPLICATIONS:  '.nreps) 
call  flout ('  MINIMUM  SOIL  CONSTANT:  ', slcmin,'  ') 

call  flout ('  MAXIMUM  SOIL  CONSTANT:  ', slcmax, '  ') 

call  flout ('  MINIMUM  ALTITUDE:  ', altmin,'  ft  ') 

call  flout ('  MAXIMUM  ALTITUDE:  ', altmax, '  ft  ') 

call  flout  ('  MINIMUM  WIND  SPEED:  ',wspmin,'  mph  ') 

call  flout  ('  MAXIMUM  WIND  SPEED:  ',wspmax, '  mph  ') 

call  flout ('  MINIMUM  WIND  DIRECTION:  ', wdrmin,'  degrees') 
call  flout  ('  MAXIMUM  WIND  DIRECTION:  ',  wdrmax, '  degrees') 


' , slcmin, ' 
' , slcmax, ' 
' , altmin, ' 
' ,  altmax, ' 
' ,  wspmin, ' 
' , wspmax, ' 


' ,  wdrmin, '  degrees' ) 
' ,  wdrmax, '  degrees ' ) 


write (3, ' (// 32x,  a/) ' ) 
write (3, ' (24x, a, ilO) ' ] 

write (3,' (//36x, a/)') 
write (3, ' (24x, a, a) ' ) 
write (3, ' (24x, a, a) ' ) 
if (minvrt .eq. 0)  then 
write (3, '  ?24x,  a) ' ) 
else 

write (3, ' (24x, a) ' ) 
endif 

call  inout ( '  NU 

call  flout ('  HE 

call  flout  ( '  B. 

call  inout ('  NUMBEF 
call  flout  ('  MINIMUM 


' **RANDOMIZATION**' 

'  MONTE  CARLO  SEED:  ', idint (rndsav) 

' **DONOR**' 

'  DESCRIPTION:  ' , dname 

'LETHALITY  DATA  FILE:  ',dfile 

'  ORIENTATION:  Erect' 

'  ORIENTATION:  Inverted' 


call  inout ('  NUMBER  OF  UNITS:  ',nunits) 

call  flout ('  HEIGHT  OF  STACK:  ',htstft,'  ft 

call  flout ('  BASE  ELEVATION:  ',htbaft,'  ft 

call  inout ('  NUMBER  OF  FRAGMENTS:  ',nfrags) 
call  flout  ('  MINIMUM  FRAGMENT  MASS:  ', fmgrmn, '  gra 
call  flout ('MINIMUM  ELEVATION  ANGLE:  ', ezndeg, '  deg 

write (3, ' (//36x,  a/) ' )  '**RANGE**' 

call  flout ('  MINIMUM  RANGE:  ', rngmin,'  ft 

call  flout ('  MAXIMUM  RANGE:  ', rngmax, '  ft 

call  inout ('  NUMBER  OF  SEGMENTS:  ',nrsegs) 

write (3, ' (//34x, a/) ' )  ' **ACCEPTOR**' 

write(3,' (20x,a,a)')  '  DESCRIPTION: 

write (3,' (20x, a, a) ' )  'VULNERABILITY  DATA  FILE: 

call  flout ('  HEIGHT  OF  STACK:  'htacft,'  ft 

call  flout ('  WIDTH  OF  STACK:  ',wdacft,'  ft 

call  flout ('  DEPTH  OF  STACK:  ',dpacft,'  ft 

close (3) 

return 


Trains  ' ) 
legrees' ) 


If  ********************************************************************** 


f ,  aname 
',afile 

■) 

') 


oooooio 

0000020 

0000030 

0000040 

0000050 

0000060 

0000070 

0000080 

0000090 

0000100 

0000110 

0000120 

0000130 

0000140 

0000150 

0000160 

0000170 

0000180 

0000190 

0000200 

0000210 

0000220 

0000230 

0000240 

0000250 

0000260 

0000270 

0000280 

0000290 

0000300 

0000310 

0000320 

0000330 

0000340 

0000350 

0000360 

0000370 

0000380 

0000390 

0000400 

0000410 

0000420 

0000430 

0000440 

0000450 

0000460 

0000470 

0000480 

0000490 

0000500 

0000510 

0000520 

0000530 

0000540 

0000550 

0000560 

0000570 

0000580 

0000590 

0000600 

0000610 

0000620 

0000630 

0000640 

0000650 

0000660 

0000670 

0000680 

0000690 

0000700 

0000710 

O000720 

0000730 

0000740 

0000750 

0000760 

0000770 

0000780 


q** ****************************  *****************************************  fooooio 


subroutine  flout (alabel, f lvar, aunit ) 

F000030 

c 

F000040 

c 

outputs  a  floating  point  variable 

F000050 

c 

with  an  alphanumeric  label  and  units 

F000060 

c 

F000070 

c 

type  declarations 

F000080 

character  alabel*25, aunit*8 

F000090 

c 

F000100 

aflvar=abs (f lvar) 

F000110 

if (aflvar.lt .1.0)  then 

F000120 

if (f lvar . It . 0 . 0 )  then 

F000130 

write (3, 9  (20x, 2a, f2 . 1, a) '  )  alabel, '  — 0' , af lvar, aunit 

F000140 

else 

F000150 

write (3, '  (2 Ox, 2a, f2 . 1, a) ' )  alabel, '  0r , f lvar, aunit 

F000160 

endif 

F000170 

else 

F000180 

write (3,' (20x, a, f 6 . 1, a) ' )  alabel, f lvar, aunit 

F000190 

endif 

F000200 

c 

F000210 

return 

F000220 

c 

F000230 

end 

F000240 

c 

F000250 

q** ********************************************************** ***********  F0002 60 

0** ************************************** *******************************  1000010 

c 

1000020 

subroutine  inout (alabel, invar) 

1000030 

c 

1000040 

c 

outputs  an  integer  variable  with  an  alphanumeric  label 

1000050 

c 

1000060 

c 

type  declaration 

1000070 

character  alabel*25 

1000080 

c 

1000090 

write (3,' (20x,a,i4)')  alabel, invar 

1000100 

c 

1000110 

return 

1000120 

c 

1000130 

end 

1000140 

c 

1000150 

0***********************************************************************  1000160 

0-** ********************************************************** ***********  A000010 

c 

A000020 

subroutine  arini (nrange, v, vd, vb, vm, value) 

A000030 

c 

A000040 

c 

initializes  variable  arrays 

A000050 

c 

A000060 

c 

local  arrays 

A000070 

dimension  v(97) , vd(97) , vb (97) , vm(97) 

A000080 

c 

A000090 

do  100  irange=l, nrange 

A000100 

v  (irange) = value 

A000110 

vd (i range) = value 

A000120 

vb (irange) =value 

A000130 

vm  (irange) rvalue 

A000140 

100 

continue 

A000150 

c 

A000160 

return 

A000170 

c 

A000180 

end 

A000190 

c  A000200 

q***********************************************************************  A000210 
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c**** ************************************************************  *******  E000010 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 


subroutine  ertim(ihbgn, tsbgn,  f rcmp, timela, uniela,  timrem, unirem) 

determines  the  elapsed  time  and  estimates  the  remaining  time 

type  declaration 
character*8  uniela, unirem 

obtain  the  current  time 
call  gettim(ihnow, imnow, isnow,  ilnow) 

compute  the  elapsed  time 
if ( ihnow . It . ihbgn )  ihnow=ihnow+2  4 
tsnow=3600 . 0*float (ihnow) +60 . 0*float (imnow) 

+  +float (isnow) +float (ilnow) /100 . 0 
timela=tsnow-tsbgn 

estimate  the  remaining  time 
t imrem=  ( 1 . 0 /f rcmp-1 . 0 ) *timela 

convert  the  time  units 
call  ctime (timela, uniela) 
call  ctime (timrem, unirem) 

return 

end 
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c***********************************************************************  E000300 


c** **************************************************************  ******* 


c 

c 

c 

c 

c 

c 


c 

c 


subroutine  ctime (time, units) 

converts  time  units 

type  declaration 
character*8  units 

if (time .gt . 60.0)  then 
time=time/60.0 
if (time .gt .60.0)  then 
time-time/ 60 . 0 
units='  hours  ' 
else 

units='  minutes' 
endif 
else 

units='  seconds' 
endif 

return 

end 


c 

c** ************************************************************ ********* 
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c: 

Ar*  **************************************  ******************************* 

•  R000010 

c 

R000020 

subroutine  rpini (rndnxt,  slcnst, asldeg, alttud, wspeed, wdrdegr 

R000030 

+  wdrrad,  jxrng,  jrkmax) 

R000040 

c 

R000050 

c 

initializes  replication  loop 

R000060 

c 

R000070 

c 

common  blocks 

R000080 

common/cnstm/  qpi,  hpi,  pi,  tpi,  raddeg,  degrad,  smlrad,  smlnum,  bignum 

R000090 

common/error/  meror, irerr, if err,merout 

R000100 

common/prmtr/  slcmin,  slcmax,  altmin,  altmax. 

R000110 

+  wspmin ,  wspmax ,  wdrmin ,  wdrmax 

R000120 

common/ rangs/  nrsegs,nrsegp 

R000130 

common/pmiss/  pn (97) ,pnd(97) ,pnb (97) ,pnm(97) 

R000140 

c 

R000150 

c 

type  declaration 

R000160 

double  precision  rndnxt 

R000170 

c 

R000180 

c 

report  call 

R000190 

if (merout .gt . 0)  write (2, '  (2x, a) ' )  'rpini  called' 

R000200 

c 

R000210 

c 

randomize  parameters 

R000220 

slcnst=rnval (slcmin, slcmax, rndnxt) 

R000230 

asldeg=10 .8*slcnst**0.38 

R000240 

alttud^rnval (altmin, altmax, rndnxt) 

R000250 

wspeed^rnval (wspmin, wspmax, rndnxt) 

R000260 

wdrdeg=rnval  (wdrmin,  wdrmax,  rndnxt) 

R000270 

wdrrad=degrad*wdrdeg 

R000280 

c 

R000290 

c 

set  integration  parameters 

R000300 

if ( (wspeed.eq. 0.0) .or. (wdrdeg.eq. 0.0))  then 

R0  00310 

jxrng=0 

R000320 

jrkmax=4 

R000330 

else 

R000340 

jxrng=l 

R000350 

jrkmax=6 

R000360 

endif 

R000370 

c 

R000380 

c 

initialize  the  miss  probabilities 

R000390 

call  arini (nrsegp, pn, pnd, pnb, pnm, 1.0) 

R000400 

c 

R000410 

c 

report  return 

R000420 

if (merout .gt . 0)  write (2, r  (2x, a)')  'returning  from  rpini' 

R000430 

c 

R000440 

return 

R000450 

c 

R000460 

end 

R000470 

c 

R000480 

c**** ************************************ ************************ ******* 

R000490 
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****** ****** ****** ****************************************** ***********  F000010 


subroutine  fnm(irep,  ifrag, rndnxt,  irico, lrange,htfrft,mfrhtl, 
+  vf rag, seldeg, selrad, aeldeg, aelrad,  cdl,  cd2, cd3, 

+  deldst) 

initializes  fragment  loop 

common  blocks 


common/ cnstu/  fgrtlb, fgmtlb, fgrtgm, ffttm, f fttmm, fmtft, f i2tc2 

common/error/  meror, irerr, iferr,merout 

common/ zones/  elzdeg, azsdeg,  azsrad 

common/donor/  nunits, htstft , htbaf t 

common/f ragd/  famil(900) ,vfifs (900) ,farma (900) , 

+  fatn (900) ,ezldeg (900) 

common/acptr/  htacft , wdacft,  dpacf t , htacfm 
common/rkvar/  rkv(7) 

type  declaration 
double  precision  rndnxt 

set  the  error  output  flag 
merout=0 

if (meror .gt . 0)  then 

if (irep.gt. irerr)  then 
if  (ifrag. gt.iferr)  then 
merout=l 
endif 
endif 
endif 

report  call 
if (merout . gt . 0 ) 

+  write (2, ' (2x,a, i3) ’ )  ' frini  called  at  fragment  ifrag 

initialize  the  ricochet  count  and  last  range  segment  index 
irico=0 
lrange=0 

randomize  the  initial  fragment  height 
htf rf t=rnval (htbaf t, htstft, rndnxt) 

set  switch  1  for  fragment  height 
if (htf rft .ge .htacft)  then 
mfrhtl=l 
else 

mf rhtl=0 
endif 

randomize  the  initial  fragment  velocity 
rndlog=sqrt (-2 . *log (rndom( rndnxt) ) ) 
rndnum=tpi*rndom ( rndnxt ) 
rndcos=rndlog*cos (rndnum) 
rndsin=rndlog* sin (rndnum) 

vf rag=vf if s (ifrag) +vf ifs (ifrag) *0 . 035*rndcos 

randomize  the  initial  fragment  elevation  angle 
seldeg= (ezldeg (ifrag) +elzdeg*rndom( rndnxt) ) 
if (seldeg. gt . 89 . 99)  then 
seldeg=89. 99 
else 

if ( (seldeg. It. 1.0e-2) .and. (seldeg. ge . 0 . 0) )  then 
seldeg=l . 0e-2 
else 

if ( ( seldeg . It . 0 . 0 ) . and . ( seldeg . gt . -1 . 0e-2 ) )  seldeg=-l . 0e-2 
endif 
endif 

selrad=degrad* seldeg 
aeldeg=abs (seldeg) 
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c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c**** 


aelrad=abs  (selrad) 

randomize  drag  coefficient  parameters 
cdmax=l . 75*farma (if rag) -1.27 
cdmin=0 . 66*farma (if rag) -0.26 
cd0=rnval (cdmin, cdmax, rndnxt ) 
cdl=cd0+0 .2 
cd2=cd0+0 . 65 
cd3=cd0+0 . 5 

initialize  the  integration  step 
call  dstep (htfrft, seldeg, selrad,  deldst) 

initialize  the  integration  variables 
rkv (1) =0 . 0 

rkv  (2)  ==vf  rag*cos  (selrad) 
rkv (3) =0.0 

rkv (4) =vf rag*sin (selrad) 

rkv(5)=htfrft 

rkv (6) =0.0 

rkv (7) =0.0 

report  return 

if (merout .gt . 0)  write (2,'  (2x,a)')  'returning  from  frini' 

return 

end 
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************************************************************ *******  pQ  01 020 


c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 


***********************************************************************  DO 000 10 


subroutine  dstep (htfrft, seldeg, selradf  deldst) 
determiness  integration  distance  step 
common  blocks 

common/error/  meror, irerr, if err, merout 

common/range/  rngmax, rngseg, rngsgh, rngscl,  rngsix,  rngmin 
report  call 

if (merout .gt . 0)  write (2, ' (4x, a) ' )  'dstep  called' 

if (seldeg.lt .0.0)  then 

dtognd=htf rft/sin (-selrad) 
if (dtognd . gt . rngsix)  then 
deldst=rngsgh 
else 

deldst =aint (dtognd/ 1.2) 
endif 
else 

if (seldeg. gt .70 .0)  then 
deldst=rngsgh 
else 

deldst=rngscl*aint { (-0. 02492*seldeg+2 .20134)  *seldeg+18 . 8306) 
endif 
endif 

report  return 

if (merout .gt .0)  write (2, ' (4x, a) ' )  ' returning  from  dstep' 

return 

end 


D000020 

D000030 

D000040 

D000050 

D000060 

D000070 

D000080 

D000090 

D000100 

D000110 

D000120 

D000130 

D000140 

D000150 

D000160 

D000170 

D000180 

D000190 

D000200 

D000210 

D000220 

D000230 

D000240 

D000250 

D000260 

D000270 

D000280 

D000290 

D000300 

D000310 

D000320 

D000330 

D000340 

D000350 


***********************************************************************  DO 003 60 


56 


oo  ooooo  oooooo  oooo  oooooooo  oooooo 


c 

c 

c 


*********************************************************************** 
subroutine  tstep (deldst, irico, vfrag, aelrad,  deltim) 
determiness  integration  time  step 
common  blocks 

common/cnstm/  qpi,hpi,pi,tpi,  raddeg, degrad,  smlrad,  smlnum, bignum 
common/error/  meror, irerr, iferr,merout 
common/acptr/  htacft, wdacft , dpacf t , htacfm 
common / rk va r /  rkv ( 7 ) 

report  call 

if (merout .gt .0)  write (2, ' (2x, a) ' )  'tstep  called' 

compute  the  standard  time  step 
deltim=deldst /vfrag 

compute  a  special  time  step  if  approaching 
the  top  of  thr  hazard  volume  or  the  ground 
if (rkv(4) .It. 0.0)  then 

elerad=max (aelrad, smlrad) 
deldsp=deldst+htacfm 

if  the  fragment  is  below  the  top  of  the  hazard  volume 
if (rkv (5) .le. htacft)  then 

confute  a  time  step  approaching  the  ground 
dtognd=rkv ( 5 ) /sin (elerad) 

if (dtognd.lt .deldsp)  deltim= (dtognd-htacfm) /vfrag 
if (dtognd.lt .htacft)  deltim= (dtognd-0 . 3) /vfrag 
if (dtognd.lt .0.4  )  deltim=dtognd/vf rag+2 . Oe-6 

else 

if  the  fragment  is  low  enough 
if (rkv(5) .It. (deldst+12 . 0) )  then 

compute  a  time  step  approaching 
the  top  of  thr  hazard  volume 
dtotop= (rkv (5) -htacft) /sin (elerad) 
if (dtotop. It .deldsp)  deltim= (dtotop-htacfm) /vfrag 
if (dtotop.lt .htacft)  deltim= (dtotop-0.3) /vfrag 
if (dtotop . It . 0 . 4  )  deltim=dtotop/vf rag+2 . Oe-6 

endif 

endif 

endif 

limit  the  time  step 

if ( (vfrag. It. 100.0) .and. (rkv(4) .gt.-30.0) .and. (deltim. gt . 0 . 1) ) 

+  deltim=0 . 1 

if ( (irico. gt .1) .and. (vfrag.lt . 150 .0) .and. (deltim. gt . 0 . 08) ) 

+  deltim=0.08 

report  return 

if  (merout. gt.0)  write (2, ' (2x, a) ' )  ' returning  from  tstep' 

return 

end 
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c***********************************************************************  R0OOO1O 


c  R000020 

subroutine  rkint  (wspeed,  wdrdeg,  wdrrad,  jxrng,  jrkmax,  cdlf  cd2,  cd3,  R000030 

+  deltim,  alttud, if rag)  R000040 

c  R000050 

c  performs  Runge-Kutta  integration  R000060 

c  R000070 

c  common  blocks  RO 00080 

common/cnstm/  qpi, hpi,pi, tpi,  raddeg, degradf  smlrad,  smlnum,  bignum  R000090 

common/cnstp/  gg,  t glbgr , hrair 0 , racnst , vsound, hccnst , hcaml , hcam2 ,  RO  0  0 1 0 0 

+  hcam3  R000110 

common/error/  meror, irerr, iferr, merout  R000120 

common/fragd/  famil(900) ,vfifs(900) , farma{900) ,  R000130 

+  fatn(90Q) ,ezldeg(900)  R000140 

common/rkvar/  rkv(7)  R000150 


c  R000160 

c  local  arrays  R000170 

dimension  rko (7) , rks (7) , rk (4,  7)  R000180 


c  R000190 

c  report  call  R000200 

if (merout .gt . 0)  write (2, ' (2x, a) ' )  ' rkint  called'  R000210 

c  R000220 

c  begin  Runge-Kutta  outer  loop  R000230 

do  200  irk=l,4  R000240 


c  R000250 

if (irk . It . 3)  then  R000260 

rkcnst=0 .5  R000270 

else  R000280 

rkcnst=l.  R000290 

endif  R000300 


c 

c  determine  relative  wind 

if ( wspeed.  le. 0.0)  then 

vfrel=sqrt (rkv(2) *rkv(2) +rkv(4) *rkv(4) ) 
aayrad=atan(rkv(4) /rkv(2) ) 
aaxrad=hpi-aayrad 
else 

if ( wdrdeg. eq. 0.0)  then 
vxf rel=rkv (2 ) -wspeed 
vyfrel=rkv(4) 

vfrel=sqrt (vxfrel*vxf rel+vyfrel*vyf rel) 
aayrad=atan(vyfrel/abs (vxf rel) ) 
aaxrad=atan (vxfrel/abs (vyf rel) ) 
else 

vxf rel=rkv (2) -wspeed*cos (wdrrad) 
vyfrel=rkv(4) 

vzf rel=wspeed*sin (wdrrad) -rkv<6) 

vf rel=sqrt (vxfrel*vxf rel+vyf rel*vyf rel+vzf rel*vzf rel) 
aayrad=atan(vyfrel/sqrt (vxf rel*vxf rel+vzf rel*vzf rel) ) 
aaxrad=atan (vxf rel/sqrt ( vyf rel*vyf rel+vzf rel*vzf rel) ) 
aazrad=atan (vzfrel/sqrt (vxf rel*vxf rel+vyf rel*vyf rel) ) 
endif 
endif 
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c  determine  air  density,  Mach  number  and  drag  coefficient 

hdelt=-( ( rk v ( 5 ) talttud) +rkv(4) *0 .5*deltim) 
hrhair=hrairO*expon (hdelt/racnst) 
hcam=vfrel/ (vsound*expon (hdelt/hccnst ) ) 
if (hcam.ge.hcaml)  then 
if (hcam.ge.hcam2)  then 
if (hcam.ge .hcam3)  then 
cd=cd3 
else 

cd=cd2-0 . 15/  (hcam3-hcam2)  *  (hcam-hcam2) 
endif 
else 

cd=cd2-0 .45/  (hcam2-hcaml)  *  (hcam2-hcam) 
endif 
else 

cd=cdl-0 .2/ (hcaml-0 . 1) * (hcaml-hcam) 
endif 
c 

c  compute  velocities  and  accelerations 

acdrag=hrhair*vf rel*vf rel*cd*f amil (ifrag) /144 . 
rko (2) =-acdrag*sin (aaxrad) 
rko (3) =rkv (2) 

rko (4) =-acdrag*sin (aayrad) -gg 

rko (5)=rkv(4) 

if ( jxrng . ne . 0 )  then 

rko (6) =acdrag*sin (aazrad) 
rko (7)=rkv(6) 
endif 
c 

c  begin  Runge-Kutta  inner  loop 

do  100  jrk=2, jrkmax, 2 

if ( .not . (irk.gt . 1) )  rks ( jrk) =rkv ( jrk) 
rk (irk, jrk) =rko ( jrk) *deltim 
if (irk. eq. 4)  then 
rkv(jrk)=rks ( jrk)+ 

+  (rk (1, jrk) +2*rk (2, jrk) +2*rk (3, jrk) +rk (4, jrk) ) /6 . 0 
rk v ( j  rk+ 1 ) =rkv ( j  r k+ 1 )  + 

+  deltim* (rks ( jrk) + (rk (1, jrk)+rk(2, jrk)+rk<3, jrk) )/ 6.0) 
else 

rkv (jrk) =rks (jrk) +rk (irk, jrk) *rkcnst 
endif 

100  continue 
c 

200  continue 
c 

c  report  return 

if (merout.gt .0)  write (2, ' (2x, a) f )  ' returning  from  rkintf 
c 

return 

c 

end 

c 
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*********************************************************************** 
subroutine  f rvar ( vf rag,  selrad, seldeg, aelrad, aeldegr  range, irange) 
determines  the  fragment  variable  values 
common  blocks 

common /cnstm/  qpi,  hpi,pi,  tpi,  raddeg, degrad,  smlrad,  smlnum, bignum 
common/error/  meror, irerr, if err , merout 
common/rangs/  nrsegs, nrsegp 

common/range/  rngmax,  rngseg, rngsgh, rngscl, rngsix, rngmin 
common/ rkvar/  rkv(7) 

report  call 

if (merout .gt .0)  write (2, ' (2x, a) '  )  'frvar  called' 
compute  the  fragment  velocity 

vfrag=sqrt (rkv(2) *rkv(2) +rkv(4) *rkv(4) +rkv(6) *rkv(6) ) 
compute  the  elevation  angle 

selrad=atan (rkv (4) /sqrt (rkv(2) *rkv (2) +rkv ( 6) *rkv(6) ) ) 

seldeg=raddeg*selrad 

aelrad=abs (selrad) 

aeldeg=abs (seldeg) 

compute  the  range  and  segment  index 
range=sqrt (rkv (3) *rkv(3) +rkv(7) *rkv(7) ) 
irange=int (range/rngseg) +1 
call  irlim (irange, nrsegs, nrsegp) 

report  return 

if (merout .gt . 0)  write (2, ' (2x, a) ' )  'returning  from  frvar' 

return 

end 

*r*******************************************************  ********  ******* 


F000010 

F000020 

F000030 

F000040 

F000050 

F000060 

F000070 

F000080 

F000090 

F000100 

F000110 

F000120 

F000130 

F000140 

F000150 

F000160 

F000170 

F000180 

F000190 

F000200 

F000210 

F000220 

F000230 

F000240 

F000250 

F000260 

F000270 

F000280 

F000290 

F000300 

F000310 

F000320 

F000330 

F000340 

F000350 

F000360 

F000370 

F000380 


***********************************************************************  1000010 

1000020 

subroutine  i r 1 im ( i range ,  nrsegs , nrsegp ) 


1000030 
1000040 
1000050 
1000060 
1000070 
1000080 
1000090 
1000100 
1000110 
1000120 
1000130 

***********************************************************************  jo 00140 


limits  the  value  of  irange 

if (irange.lt. 1)  irange=l 

if ( irange. gt .nrsegs)  irange=nrsegp 

return 

end 
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no  o  oooo  oooooo 


c 

c 

c 


*********************************************************************** 
subroutine  hazch  (if rag,  irange, mfrhtl,  mfrht2f  Irange,  aelrad, vfrag) 
directs  hazard  computations  if  they  are  required 
common  blocks 

common/error/  meror, irerr, iferr, merout 
common/ rangs /  nrsegs , nrsegp 

common/ range/  rngmax,  rngseg,  rngsgh,  rngscl, rngsix, rngmin 
common/acptr/  htacft, wdacft, dpacftf htacfm 
common/rkvar/  rkv(7) 

report  call 

if (merout .gt .0)  write (2, ' (2x, a) f )  'hazch  called' 

set  switch  2  for  fragment  height 
if (rkv(5) .gt. htacft)  then 
mf  rht  2=1 
else 

mf  rht 2=0 
endif 

if (mfrht2 .eq. 1)  then 
if (mfrhtl .ne. 1)  then 
mf rht 1=1 

delrng= (rkv (5) -htacft) /tan (aelrad) 
irange= 

+  int ( (sqrt (rkv(3) *rkv(3) +rkv(7) *rkv(7) )-delrng) /rngseg) +1 
call  irlim(irange, nrsegs, nrsegp) 
call  hazrd(ifrag, irange, Irange,  aelrad,  vfrag) 
endif 
else 

call  hazrd (if rag, irange, Irange, aelrad, vfrag) 
endif 

report  return 

if  (merout .gt . 0)  write (2,'  (2x,a)r)  'returning  from  hazch' 

return 

end 


H000010 

H000020 

H000030 

H000040 

H000050 

H000060 

H000070 

H000080 

H000090 

H000100 

H000110 

H000120 

H000130 

H000140 

H000150 

H000160 

H000170 

H000180 

H000190 

H000200 

H000210 

H000220 

H000230 

H000240 

H000250 

H000260 

H000270 

H000280 

H000290 

H000300 

H000310 

H000320 

H000330 

H000340 

H000350 

H000360 

H000370 

H000380 

H000390 

H000400 

H000410 

H000420 

H000430 


q**** ****************************************************** *************  HO 00440 
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***********************************************************************  HO o 0010 


c 

c 


subroutine  hazrd(if rag, irange, Irange,  aelrad, vfrag) 
performs  hazard  computations 
common  blocks 

common/cnstm/  qpi,hpi,pi,  tpi,  raddeg,  degrad,  smlrad,  smlnum, bignum 

common/ error/  meror, irerr, iferr,merout 

common/ zones/  elzdeg, azsdeg, azsrad 

common/f ragd/  famil (900) , vfifs (900) ,farma (900) , 

+  fatn (900) ,ezldeg(900) 

common/ range/  rngmax,  rngseg, rngsgh, rngscl, rngsix, rngmin 
common/arang/  rangem(97) ,ahfrn(97) ,ahtop(97) 
common/ acpt r/  htacf t , wdacf t , dpacf t , ht acfm 
common/ acf lg/  nweps , inacwh,  macrm,  macsh ,  mstor 
common/thorw/  athw, bthw, cthw, dthw, vrwms ,  tew 
common/thorr/  athr , bthr , cthr , dthr, vrrms , ter 
common/thras/  athasv,bthasv, cthasv, tcasv, 

+  athasm, bthasm, cthasm, tcasm 

common/vulns/  amf rn, amtop/axf  rn, axtop, apf rn, aptop 
common/ rkvar/  rkv (7) 

common/pmiss/  pn (97) ,pnd(97) ,pnb (97) ,pnm(97) 
common/denst/  fd(97) , fdd<97) , fdb (97) ,  fdm(97) 
common/phitt/  ph(97) ,phd(97) ,phb(97) ,phm(97) 

report  call 

if (merout . gt . 0 )  write (2, 9  (4x, a) ' )  'hazrd  called' 

initialize  parameters  if  range  segment  index  has  changed 
if ( i range. ne. Irange)  then 

call  wiini  (itrpt , pnO ,  pnOd, pnOb,  pnOm,  pnr , pnrd, pnrb,  pnrm, 

+  fdO,  fdOd,  fdOb,  fdOm,  fdr,  fdrd,  fdrb,  fdrm) 

endif 

increment  within  segment  counter 
it  rpm=  itrpt 
itrpt-itrpt+1 

compute  the  sine  and  cosine  of  the  elevation  angle 
sinele=sin (aelrad) 
cosele=cos (aelrad) 

compute  the  hazard  volume  presented  area 
apvft2=ahfrn (irange) *cosele+aht op (irange) *sinele 

compute  the  total  presented  area  of  the  acceptor 
apaft2=htacft*cosele 

if (rkv (4) .It. 0.0)  apaft2=apaft2+dpacft*sinele 
apaft2=apaft2*wdacft 

check  and  adjust  fragment  parameters  for  shroud  penetration 
call  shpen ( if rag, vfrag,  vf haz , vf cdx, vf cdp, vf cbx, vf cbp) 

accumulate  the  unconditional  miss  probability 
call  nohit (if rag, apaft2, apvft2, irange, itrpm, itrpt, pnr, pn0,pn) 

accumulate  the  unconditional  fragment  density 
call  dnsty (if rag, apvft2, irange, itrpm,  itrpt, fdr, fdO, fd) 


H000020 

H000030 

H000040 

H000050 

H000060 

H000070 

H000080 

H000090 

H000100 

H000110 

H000120 

H000130 

H000140 

H000150 

H000160 

H000170 

H000180 

H000190 

H000200 

H000210 

H000220 

H000230 

H000240 

H000250 

H000260 

H000270 

H000280 

H000290 

H000300 

H000310 

H000320 

H000330 

H000340 

H000350 

H000360 

H000370 

H000380 

H000390 

H000400 

H000410 

H000420 

H000430 

H000440 

H000450 

H000460 

H000470 

H000480 

H000490 

H000500 

H000510 

H000520 

H000530 

H000540 

H000550 

H000560 

H000570 

H000580 

H000590 
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c  determine  the  detonation  vulnerable  area 

avdf t2=0 . 0 
if (macwh . ne . 0 )  then 

call  vulnd(vfcdx, vfhaz,axfrn, axtop, cosele, sinele,  avdft2) 
endif 

if  (macrm.ne.  0)  then 

call  vulnd(vfcdp, vfhazr apfrn, aptop, cosele, sinele,  avdft2) 
endif 
c 

c  accumulate  the  detonation  miss  probability  and  fragment  density 
call  nohit  (if rag,  avdft2,apvft2,  irange,  itrpm,  itrpt,  pnrd,pnOd,pnd) 
if (avdft2 .ne . 0 . 0)  then 

call  dnsty (if rag,  apvft2, irange , itrpm, itrpt,  fdrd,  fdOd, fdd) 
endif 
c 

c  determine  the  burn  vulnerable  area 

avbft2=0 . 0 
if  (macwh . ne . 0 )  then 

call  vulnb (vf cbx, cthw,  dthw, vrwms, vf haz, axf rn, axtop,  cosele, 

+  sinele,  avbft2) 

endif 

if (macrm.ne .0)  then 

call  vulnb ( vf cbp,  cthr, dthrf  vrrms , vf haz, apf rn, aptop, cosele, 

+  sinele,  avbft2) 

endif 
c 

c  accumulate  the  burn  miss  probability  and  fragment  density 

call  nohit  (if rag, avbft 2,  apvft2,  irange,  itrpm,  itrpt,pnrb,pnOb,pnb) 
if (avbft2 .ne . 0 . 0)  then 

call  dnsty (if rag, apvft 2, irange, itrpm, itrpt , fdrb, fdOb, fdb) 
endif 


c 

c 

c 

c 

c 


c 

c 

c 

c 


determine  the  mechanical  damage  vulnerable  area 
call  vulnm(ifrag, vfhaz, cosele, sinele,  avmft2) 

accumulate  the  mechanical  damage 
miss  probability  and  fragment  density 

call  nohit  (if rag,  avmft2,  apvft2,  irange,  itrpm,  itrpt, pnrm,  pn0m,pnm) 
if  (avmft2 .ne . 0 . 0)  then 

call  dnsty (ifrag, apvft2, irange, itrpm, itrpt, fdrm,  fdOm, fdm) 
endif 

lrange=i range 
report  return 

if (merout .gt . 0)  write (2, f (4x,a)')  'returning  from  hazrd' 
return 


c 

end 

c 

c******** ****************************** ****** ****** ********** **** ******* 
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H000620 
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H000650 
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H000960 

H000970 

H000980 

H000990 
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H001020 

H001030 

H001040 

H001050 

H001060 

H001070 

H001080 

H001090 

H001100 

H001110 

H001120 
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Q*********************************************************************** 

c 

subroutine  wiini (it rpt , pnO , pnOdf pnOb, pnOm, pnr , pnrd, pnrb, pnrm, 

+  f dO f  f  dOd, f dOb, f dOm,  f dr ,  f  drd,  f drb, f drm) 

c 

c  initializes  within  segment  parameters 

c 

c  common  block 

common/error/  meror, irerr, if err, merout 
c 

c  report  call 

if (merout. gt.O)  write (2,  '  <6x, a) '  )  'wiini  called' 

initialize  the  within  segment  point  index 
itrpt=0 

initialize  the  miss  probability  parameters 
pnO=l . 0 
pnr=0 . 0 
pnOd=l . 0 
pnrd=0 . 0 
pnOb=l . 0 
pnrb=0 . 0 
pnOm=l . 0 
pnrm=0 . 0 

initialize  the  fragment  density  parameters 
f  dr=0 . 0 
fd0=0 . 0 
fdrd=0 . 0 
fdOd-O . 0 
f drb=0 . 0 
fd0b=0 . 0 
fdrm=0 . 0 
fd0m=0 . 0 

report  return 

if (merout. gt.O)  write (2 , ' (6x,a)')  'returning  from  wiini' 
return 
end 


**** ************ ★*★*★****•*★***★* **************************** *********** 


wooooio 
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***********************************************************************  so 000 10 

^  S000020 

subroutine  shpen (ifrag, vf rag,  vfhaz, vfcdx, vfcdp, vfcbx, vfcbp)  S000030 


subroutine  shpen (if rag,vfrag,  vfhaz, vfcdx, vfcdp, vfcbx, vfcbp)  S000030 

w  ,  ^  ^  S000040 

checks  and  adjusts  fragment  parameters  for  shroud  penetration  S000050 

_  ,  S000060 

common  blocks  S000070 

common/cnstm/  gpi,  hpi,j?i,  tpi,  raddeg,  degrade  smlrad,  smlnum, bignum  S000080 
common/cnstu/  fgrtlb,  fgmtlb,  fgrtgm,  ffttm,  f  fttmm,  fmtft,  f i2tc2  SO 000 90 

common/error/  meror, irerr, iferr,merout  S000100 

common / f r agd/  famil(900) ,vfifs (900) , farma(900) ,  S000110 

+  fatn(900) ,ezldeg(900)  S000120 

common / f rage /  facm2  (900)  , fdmm(900) , fmgm(900)  S000130 

common /fragx/  vcrdx(900) , vcrbx(900)  S000140 

common/fragp/  verdp (900) f verbp (900)  S000150 

common/acf lg/  nweps,macwh,macrm,macsh,mstor  SO 001 60 

common/geown/  whim, whdin, whhin, whhem, whhmm, exlin  S000170 

common/ j rdtx/  abjrx, cjrx  SO 00 180 

common/thorw/  atnw, btnw? cthw, dthw, vrwms ,  t  cw  S000190 

common/georm/  rmlin,  rmdin,  rmhin,  rmhem,  rmhmm, prlin  S000200 

common/ jrdtp/  abjrp,cjrp  S000210 

common/thorr/  athr, btnr, cthr, dthr, vrrms, ter  S000220 

common/geoas/  ashem  S000230 

common/thras/  athasv,bthasv, cthasv, tcasv,  S000240 

+  athasm,bthasm,  cthasm,  tcasm  S000250 

S000260 

report  call  S000270 

if (merout.gt .0)  write (2, ' (6x, a) ' )  'shpen  called'  S000280 

S000290 

if (macsh.eq. 0)  then  S000300 

.  ^  .  S000310 

no  shroud  is  present  S000320 

^  S000330 

assign  unadjusted  values  to  hazard  fragment  parameters  S000340 

vfhaz=vfrag  S000350 

vfcdx=vcrdx (ifrag)  S000360 

vf cdp=vcrdp (ifrag)  S  0  0  0  37  0 

vf cbx=vcrbx (ifrag)  S000380 

vfcbp=vcrbp (ifrag)  S000390 

,  S000400 

else  S000410 

^  J  .  S000420 

a  shroud  is  present  S000430 

^  ^  .  S000440 

adjust  the  fragment  velocity  and  mass  for  shroud  penetration  SO 00 450 

call  vmthr (ashem,  facm2 (ifrag) , vf rag, fmgm(ifrag) ,tcasv,  athasv,  S000460 
+  bthasv, cthasvf tcasm, athasm, bthasm, cthasm,  vfhaz,  SO 00470 

+  fmhaz)  S000480 

^  ,  S000490 

compute  the  hazard  fragment  area  in  sq  cm  S000500 

fhacm2=fi2tc2*famil (ifrag) *fgmtlb*fmhaz  S000510 

^  ^  .  .  .  S000520 

compute  the  hazard  fragment  diameter  in  mm  S000530 

fhdmm=2 . 0*sqrt  <100 . 0*fhacm2/pi)  S000540 

^  .  S000550 

compute  the  explosive  detonation  critical  velocity  SO 005 60 

vf cdx=vcr jr (ab^rx, cjrx, whhmm, fhdmm)  S000570 

S000580 

compute  the  propellant  detonation  critical  velocity  S000590 

vfcdp=vcr jr (ab jrp, c jrp, rmhmm, fhdmm)  S000600 

,  .  .  S000610 

compute  the  explosive  burning  critical  velocity  S000620 

vfcbx=vcrth (tew, athw,bthw, whhem,  f ha cm2, fmhaz)  S000630 

S000640 

compute  the  propellant  burning  critical  velocity  S000650 

vfcbp=vcrth (tcr,athr,bthr, rmhem,  fhacm2, fmhaz)  S000660 

S000670 

endif  S000680 

S000690 

.report  return  S000700 

if (merout.gt .0)  write(2, ' (6x,a) ' )  'returning  from  shpen'  S000710 

S000720 

return  S000730 

S000740 

end  S000750 

S000760 

***********************************************************************  s 00 0770 


a  shroud  is  present 


f  fmhaz) 

compute  the  hazard  fragment  area  in  sq  cm 
fhacm2=f i2tc2*famil  (ifrag)  *fgmtlb*fmhaz 

compute  the  hazard  fragment  diameter  in  mm 
fhdmm=2 . 0*sqrt (100 . 0*fhacm2/pi) 

compute  the  explosive  detonation  critical  velocity 
vf cdx=vcr jr (abjrx, cjrx, whhmm, fhdmm) 

compute  the  propellant  detonation  critical  velocity 
vfcdp=vcr  jr  (ab  jrp,  c  jrp,  rmhmm,  fhdmm) 

compute  the  explosive  burning  critical  velocity 
vfcbx=vcrth (tew, athw,bthw, whhem, fhacm2, fmhaz) 

compute  the  propellant  burning  critical  velocity 
vfcbp=vcrth (ter, athr,bthr, rmhem,  fhacm2, fmhaz) 

endif 

t report  return 

if (merout.gt .0)  write(2, '  (6x,a) ' )  'returning  from  shpen' 
return 
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***********************************************************************  V00001Q 

V000020 

subroutine  vmthr  (shhcm,  facm2,  vf  fs,  fmgm,  tcsv,  athsv,bthsv,  cthsv,  VO 00030 

+  tcsm, athsm,bthsmr cthsm,  vf,fm)  V000040 

V000050 

adjusts  the  fragment  velocity  and  mass  for  shroud  penetration  V000060 

V000070 

common  blocks  V000080 

common/cnstu/  fgrtlb,  fgmtlb,  fgrtgm,  f  fttm,  f fttmm,  fmtft,  f i2tc2  V000090 

common/error/  meror, irerr, iferr, merout  V000100 

V000110 

report  call  V000120 

if (merout .gt.O)  write (2, '  <8x, a) '  )  'vmthr  called'  V000130 

V000140 

compute  the  plug  volume  in  sq  cm  VO 00150 

ha=snhcm*  f acm2  V0  0  0 1 6  0 

V000170 

convert  the  fragment  velocity  from  ft/s  to  m/s  V000180 

vf  ms=f  f  ttm*vf  f  s  V0  0  0 1 90 

V000200 

adjust  the  fragment  velocity  in  m/s  using  THOR  (zero-obliquity)  V000210 
vf=vfms-tcsv* (ha**athsv) * (fmgm**bthsv) * (vfms**cthsv)  V000220 

if (vf.lt. 0.0)  vf-0.0  V000230 

V000240 

adjust  the  fragment  mass  in  gm  using  THOR  (zero-obliquity)  V000250 

fm=fmgm-tcsm*  (ha**athsm)  *  (fmgm**bthsm)  *  (vfms**cthsm)  V000260 

if  (fm.lt. 0.0)  fm==0 . 0  V000270 

V000280 

convert  the  hazard  fragment  velocity  from  m/s  to  ft/s  V000290 

vf=fmtft*vf  V000300 

V000310 

report  return  V0 00 320 

if (merout .gt . 0)  write(2, ' (8xfa)f)  'returning  from  vmthr'  V000330 

V000340 

return  V000350 

V000360 

end  V000370 

V000380 

***********************************************************************  V0 003 90 

***********************************************************************  V000010 

V000020 

subroutine  vulnd (vcrd, vfhaz, af rn, atop, cosele, sinele,  avdft2)  V000030 


subroutine  vulnd (vcrd, vfhaz, af rn, atop, cosele, sinele,  avdft2)  V000030 

V000040 

determines  the  detonation  vulnerable  area  V000050 

V000060 

common  blocks  V0 00 070 

common/cnstu/  fgrtlb, fgmtlb, fgrtgm, f fttm, f fttmm, fmtft, f i2tc2  V000080 

common/error/  meror, irerr, if err, merout  V000090 

V000100 

report  call  V000110 

if (merout .gt.O)  write (2, ' (6x, a) ' )  'vulnd  called'  V000120 

V000130 

convert  the  fragment  velocity  from  ft/s  to  mm/us  V000140 

vfmmus=f  fttmm*  vfhaz  V0  00 150 

V000160 

check  the  detonation  condition  V000170 

if (vfmmus .gt . vcrd)  then  V000180 

V000190 

the  fragment  can  produce  detonation  V000200 

V000210 

compute  the  maximum  obliquity  V000220 

cosobx=vcrd/ vfmmus  V0  00230 

V000240 

compute  the  vulnerable  area  for  initiation  of  detonation  V000250 

call  iniar (cosobx, cosele, sinele, af rn, atop,  avdft2)  V000260 

V000270 

endif  V0 002 80 

V000290 

report  return  V0  00300 

if (merout .gt.O)  write (2,' (6x,a)')  'returning  from  vulnd'  V000310 

V000320 

return  V000330 

V000340 

end  V000350 

V000360 

***********************************************************************  v0 00370 
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o  o 


0 ************************************** ************ *********************  VO 000 10 
c  V000020 

subroutine  vulnb  (vcrb,cth,dth, vrms, vfhaz,  af  rn,atop,cosele,sinele,  VO 00030 
+  avbf t2 )  V000040 

c  V000050 

c  determines  the  burn  vulnerable  area  V000060 

c  V000070 

c  common  blocks  V0 000 80 

common/cnstu/  fgrtlb,  fgmtlb,  fgrtgm, f fttm, ffttmm, fmtft, f i2tc2  V000090 

common/error/  meror, irerr, iferr,merout  V000100 

c  V000110 

c  report  call  V000120 

if (merout.gt.0)  write (2, ' (6x, a) ' )  'vulnb  called'  V000130 

c  V000140 

c  convert  the  fragment  velocity  from  ft/s  to  m/s  V000150 

vfms=f fttm*vfhaz  V0 001 60 

c  V000170 

c  compute  the  burn  critical  condition  using  THOR  V000180 

vright=vcrb*vfms**cth  V000190 

vleft=vfms-vrms  V000200 

V000210 

check  the  burn  condition  V000220 

if (vleft.gt.vright)  then  V0 00230 

c  V000240 

c  the  fragment  can  produce  burning  V000250 

c  V000260 

c  compute  the  maximum  obliquity  V0 00270 

cosobx= (vright/vleft) **dth  V0 002 80 

c  V000290 

c  compute  the  vulnerable  area  for  initiation  of  burning  V000300 

call  iniar (cosobx, cosele,sinele,afrn,atop,  avbft2)  V000310 

c  V000320 

endif  V000330 

c  V000340 

c  report  return  V0 00350 

if (merout . gt . 0 )  write (2, ' (6x, a) ' )  ' returning  from  vulnb'  V000360 

c  V000370 

return  V000380 

c  V000390 

end  V000400 

c  V000410 

q********** ************************************************** ***********  V000420 
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oo  o  o  no  o  o  oo  oo  no  noon  o  o  o  on  ooo  oo  o  oo  on  noon  o  no  noon  on 


**************************** ********************** ************** 
subroutine  iniar (cosobx, cosele, sinele,afrnf atop,  avft2) 
computes  the  vulnerable  area  for  initiation 
common  blocks 

common/error/  meror, irerr, iferr,merout 
common/acflg/  nweps,macwh,macrm,macsh,mstor 
common/rkvar/  rkv(7) 

report  call 

if  (merout.gt.O)  write(2, '  (8x,a) ' )  'iniar  called' 

if  (mstor .eq. 1)  then 

the  storage  arrangement  is  vertical 

compute  the  maximum  tangency 
costax=cosobx/cosele 

if  the  maximum  obliquity  exceeds  the  elevation  angle 
if (costax.lt. 1.0)  then 

compute  the  front  vulnerable  area 
avft2=avft2+cosele*sqrt (1.0-costax*costax) *afrn 

endif 

if  the  vertical  fragment  velocity  component  is  negative 
if ( rkv ( 4 ) . It . 0 . 0 )  then 

and  if  the  elevation  angle  exceeds  the 
complement  of  the  maximum  obliquity 
if (sinele.gt. cosobx)  then 

compute  the  top  vulnerable  area 
avft2=avft2+sinele*atop 

endif 

endif 

else 

the  storage  arrangement  is  horizontal 

compute  the  sine  of  the  maximum  obliquity 
sinobx=sqrt (1 . 0-cosobx*cosobx) 

compute  the  front  vulnerable  area 
avft2=avft2+sinobx*af rn 

if  the  vertical  fragment  velocity  component  is  negative 
if (rkv(4) .lt.0.0)  then 

compute  the  top  vulnerable  area 
avft2=avft2+sinobx*atop 

endif 

endif 

report  return 

if  (merout.gt.O)  write (2,' (8x, a)')  'returning  from  iniar' 
return 


*******  1000010 
1000020 
1000030 
1000040 
1000050 
1000060 
1000070 
1000080 
1000090 
1000100 
1000110 
1000120 
1000130 
1000140 
1000150 
1000160 
1000170 
1000180 
1000190 
1000200 
1000210 
1000220 
1000230 
1000240 
1000250 
1000260 
1000270 
1000280 
1000290 
1000300 
1000310 
1000320 
1000330 
1000340 
1000350 
1000360 
1000370 
1000380 
1000390 
1000400 
1000410 
1000420 
1000430 
1000440 
1000450 
1000460 
1000470 
1000480 
1000490 
1000500 
1000510 
1000520 
1000530 
1000540 
1000550 
1000560 
1000570 
1000580 
1000590 
1000600 
1000610 
1000620 
1000630 
1000640 
1000650 
1000660 
1000670 


end 


1000680 

1000690 


1000700 

r*** *******************************************************************  10007 10 
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c***** *********************************************************  ********* 
c 

subroutine  vulnm(if rag, vfhaz, cosele, sinele,  avmft2) 
c 

c  determines  mechanical  damage  vulnerable  area 

c 

c  common  blocks 

common/error/  meror , irerr , if err, merout 
common/fragk/  frkec(900) 
common/mechd/  wpke 

common/ vulns /  amf  rn, amtop, axf rn, axtop, apf rn, aptop 
common/ rkvar/  rkv(7) 
c 

c  report  call 

if (merout. gt.O)  write (2, ' (6x, a) ' )  'vulnm  called' 
c 

c  compute  the  kinetic  energy 

f rke=f rkec (if rag) *vfhaz*vfnaz 
c 

.c  check  the  mechanical  damage  condition 

if (wpke.lt.frke)  then 
c 

c  the  fragment  produces  mechanical  damage 

avmft2=amfrn*cosele 

if (rkv(4) .It. 0.0)  avmft2=avmft2+amtop*sinele 
c 

else 

c 

c  the  fragment  does  not  produces  mechanical  damage 

avmf  t2=CT.  0 
c 

endif 


V000010 

V000020 

V000030 

V000040 

V000050 

V000060 

V000070 

V000080 

V000090 

V000100 

vooono 

V000120 

V000130 

V000140 

V000150 

V000160 

V000170 

V000180 

V000190 

V000200 

V000210 

V000220 

V000230 

V000240 

V000250 

V000260 

V000270 

V000280 

V000290 

V000300 

V000310 

V000320 


c 

c  report  return 

if (merout .gt . 0)  write (2, ' (6x,a)')  'returning  from  vulnm' 

return 

c 

end 

c 

c** ********************************************************** ******** 


V000330 
V000340 
V000350 
V000360 
V000370 
V000380 
V000390 
V000400 
***  V000410 


c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 


***********************************************************************  N000010 


subroutine 

+ 


nohit (if rag, aptft2, apvft2, irange, itrpm, itrpt 
pnt) 


,pnrt,pn0t. 


N000020 

N000030 

N000040 


accumulates  miss  probability 
common  blocks 

common/cnstm/  qpi,  hpi, pi,  tpi,  raddeg,  degrad,  smlrad,  smlnum,  bignum 
common/error/  meror, irerr, if err, merout 
common/fragd/  famil (900) , vfifs (900) , farma (900) , 

+  fatn (900) ,ezldeg (900) 

local  array 
dimension  pnt (97) 

report  call 

if (merout. gt.O)  write (2, ' (6x, a) ' )  'nohit  called' 
vfarat=-fatn (if rag) *aptft2/apvft2 

pnrt= (pnrt*f loat (itrpm) +expon (vfarat) ) /float (itrpt) 
pnt (irange) =max (pnt (irange) *pnrt/pn0t, 0 .0) 
pn0t=pnrt 

report  return 

if (merout. gt.O)  write (2, ' (6x,a) ' )  'returning  from  nohit' 

return 

end 


N000050 

N000060 

N000070 

N000080 

N000090 

N000100 

N000110 

N000120 

N000130 

N000140 

N000150 

N000160 

N000170 

N000180 

N000190 

N000200 

N000210 

N000220 

N000230 

N000240 

N000250 

N000260 

N000270 

N000280 

N000290 

N000300 


c 

c** ********************************************************** ****** 


N000310 
*****  N000320 
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******************** ****** *********************************************  D000010 

D000020 

subroutine  dnsty  (ifrag,  apvft 2, irange,  itrpm,  itrpt,  fdrt,  fdOt, fdt)  D000030 
accumulates  fragment  density 
common  blocks 

common/error/  meror, irerrf if err, merout 
common/ repls /  nreps , f nreps 

common/ fragd/  famil  (900) ,  vfifs  (900) ,  farma  (900) , 

+  fatn (900) , ezldeg(900) 

local  array 
dimension  fdt (97) 

report  call 

if (merout .gt . 0)  write (2, '  (6x, a) 9  )  'dnsty  called' 


fdrt= (fdrt*float (itrpm) +f atn (ifrag) /apvft2) /float (itrpt)  D000190 

fdlt=fdrt/ fnreps  D000200 

fdt (irange) =fdt (irange) +fdlt-fd0t  D000210 

fd0t=fdlt  D000220 

D000230 

report  return  D000240 

if (merout .gt . 0)  write (2f '  (6x,  a) ' )  'returning  from  dnsty'  D000250 

D000260 

return  D000270 

D000280 

end  D000290 

D000300 


************************************************** *****************  D000310 


D000040 

D000050 

D000060 

D000070 

D000080 

D000090 

D000100 

D000110 

D000120 

D000130 

D000140 

D000150 

D000160 

D000170 

D000180 
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c** ********************************************************** ***********  rooooio 


c 

subroutine  ricoc (slcnst , vf rag, aeldeg, selrad, vrico,  ricrad) 
c 

c  determines  ricochet  conditions 

c 

c  common  block 

common/ error/  meror, irerr , if err,  merout 
c 

c  report  call 

if (merout .gt . 0)  write (2, '  (2x, a) ' )  'ricoc  called' 
c 

vratl=-0 . 01597*aeldeg*aeldeg+0 . 02156*aeldeg+0 . 96170 
vrat2=-0 . 00861*aeldeg*aeldeg+0 . 00692*aeldeg+0 . 96302 
vrat3=-0 . 00387*aeldeg*aeldeg-0 . 00414*aeldeg+0 . 95592 
vrat4=-0 . 00342*aeldeg*aeldeg-0 . 00097*aeldeg+0 . 94090 
vrat5=-0 . 00243*aeldeg*aeldeg-0 . 00520*aeldeg+0 . 93080 
vrat6=-0 . 00188*aeldeg*aeldeg-0 . 00821*aeldeg+0 . 93802 
c 

eratl=  0 . 13829*aeldeg*aeldeg-0 . 98645*aeldeg+2 . 81550 
erat2=  0 . 08549*aeldeg*aeldeg-0 .78423*aeldeg+2 . 90120 
erat3=  0 . 07515*aeldeg*aeldeg-0 .73919*aeldeg+3 . 10560 
erat4=  0 . 02142*aeldeg*aeldeg-0 .37397*aeldeg+2 .78580 
erat5=  0 . 01707*aeldeg*aeldeg-0 . 32521*aeldeg+2 . 80920 
erat6=  0 . 01369*aeldeg*aeldeg-0 .29580*aeldeg+2 . 82620 
c 

if ( (slcnst. ge. 0.25) .and. (slcnst .It . 0 . 5) )  then 
vrat=vratl+ (slcnst— 0.25) /0.25* (vrat2-vratl) 
erat=eratl+ (slcnst-0 .25) /0.25* (erat2-eratl) 
go  to  100 
endif 
c 

if ( (slcnst. ge. 0.5) .and. (slcnst. It. 1.0) )  then 
vrat=vrat2+ (slcnst-0 .5) /0 .5* (vrat3-vrat2) 
erat=erat2+ (slcnst-0 .5) /0 .5* (erat3-erat2) 
go  to  100 
endif 
c 

if ( (slcnst. ge. 1.0) .and. (slcnst. It. 2.0) )  then 
vrat=vrat3+ (slcnst-1 . 0) * (vrat4-vrat3) 
erat=erat3+ (slcnst-1 . 0) * (erat4-erat3) 
go  to  100 
endif 
c 

if ( (slcnst. ge. 2.0) .and. (slcnst. It. 3.0) )  then 
vrat=vrat4+ (slcnst-2 . 0) * (vrat5-vrat4) 
erat=erat4+ (slcnst-2 . 0) * (erat5-erat4) 
go  to  100 
endif 
c 

vrat=vrat5+ (slcnst-3 . 0) * (vrat6-vrat5) 
erat=erat5+(slcnst-3 . 0) * (erat6-erat5) 
c 

100  vrico=vrat*vf rag 
ricrad=-selrad*erat 
c 

c  report  return 

if  (merout. gt.0)  write (2, '  (2x, a) ' )  ' returning  from  ricoc' 
c 

return 


c 

end 


c 

c****** ******************************************************** ********* 


R000020 

R000030 
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R000050 

R000060 

R000070 
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R000100 

R000110 

R000120 
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R000180 

R000190 

R000200 

R000210 

R000220 

R000230 

R000240 

R000250 

R000260 

R000270 

R000280 

R000290 

R0  00300 

R000310 

R000320 

R000330 

R000340 

R000350 

R000360 

R000370 

R000380 

R000390 

R000400 

R000410 

R000420 

R000430 

R000440 

R000450 

R000460 

R000470 

R000480 

R000490 

R000500 

R000510 

R000520 

R000530 

R000540 

R000550 

R000560 

R000570 

R000580 

R000590 

R000600 

R000610 

R000620 

R000630 

R000640 
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********************************** ************ ****************** ******* 

subroutine  reini  (irico,  vrico,  ricrad,  wspeed,  wdrdeg,  vf rag,  selrad, 

+  seldeg,  aelrad,  aeldeg, deldst ,  mf  rhtl) 

reinitializes  the  velocity  integration  variables  after  ricochet 

common  blocks 

common /cnstm/  qpi,  hpi,pi,  tpi,  raddeg,  degrad,  smlrad,  smlnum,bignum 
common/error/  meror, irerr, if err, merout 
common / rk va r /  r k v ( 7 ) 

report  call 

if (merout .gt . 0)  write (2,  '  (2x, a) 9  )  f reini  called' 

increment  the  ricochet  count 
irico=irico+l 

assign  the  ricochet  values  to  the  fragment  variables 
vfrag=vrico 
selrad=ricrad 
seldeg=raddeg*selrad 
aelrad=abs (selrad) 
aeldeg=abs (seldeg) 

assign  the  fragment  variable  values  to  the  Runge-Kutta  variables 
rkv (4) =vf rag*sin (selrad) 

if ( (wdrdeg.eq. 0.0) .or. (wspeed.eq. 0.0))  then 
rkv (2) =vf rag*cos (selrad) 
rkv (6) =0.0 
else 

angrad=atan(rkv(2) /rkv (6) ) 
rkv (2) =vf rag*cos (selrad) *sin (angrad) 
rkv (6) =vf rag*cos (selrad) *cos (angrad) 
endif 

call  dstep(rkv(5) , seldeg, selrad,  deldst) 
mfrhtl=0 
report  return 

if (merout .gt . 0)  write (2, '  (2x,a) 9 )  'returning  from  reini' 

return 

end 


R000010 

R000020 

R000030 

R000040 

R000050 

R000060 

R000070 

R000080 

R000090 

R000100 

R000110 

R000120 

R000130 

R000140 

R000150 

R000160 

R000170 

R000180 

R000190 

R000200 

R000210 

R000220 

R000230 

R000240 

R000250 

R000260 

R000270 

R000280 

R000290 

R000300 

R000310 

R000320 

R000330 

R000340 

R000350 

R000360 

R000370 

R000380 

R000390 

R000400 

R000410 

R000420 

R000430 

R000440 

R000450 

R000460 

R000470 


************************************** *********************************  R0 00480 
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c**** ******************************** ******  ******************** *********  HO 0  0010 
c  H000020 


subroutine  hpacc 

H000030 

c 

H000040 

c 

accumulates  hit  probabilities 

H000050 

c 

H000060 

c 

common  blocks 

H000070 

common/error/  meror, irerr, iferr, merout 

H000080 

common/ repls/  nreps, fnreps 

H000090 

common/rangs/  nrsegsf nrsegp 

H000100 

common/phitt/  ph(97) , phd(97) , phb (97) , phm(97) 

H000110 

common/pmiss/  pn(97) , pnd(97) ,pnb(97) ,pnm(97) 

H000120 

c 

H000130 

c 

report  call 

H000140 

if (merout .gt . 0)  write (2, '  (2x, a) ' )  'hpacc  called' 

H000150 

c 

H000160 

do  100  irange=l, nrsegs 

H000170 

ph  (irange)=ph  (irange) + (1 . 0-pn  (irange) ) /fnreps 

H000180 

phd(irange)=phd(irange)+(1.0-pnd(irange) ) /fnreps 

H000190 

phb  (irange) -phb (irange)  +  (1 . 0-pnb (irange) ) /fnreps 

H000200 

phm( irange) =phm( irange) + (1 . 0-pnm( irange) ) /fnreps 

H000210 

100  continue 

H000220 

c 

H000230 

c 

report  return 

H000240 

if (merout .gt . 0)  write (2 ,f  (2x,a)')  'returning  from  hpacc' 

H000250 

c 

H000260 

return 

H000270 

c 

H000280 

end 

H000290 

H000300 

***********************************************************************  HO 003 10 
***********************************************************************  O000010 


0000020 

subroutine  outpt  0000030 

0000040 

writes  output  to  appropriate  files  0000050 

0000060 

common  blocks  0000070 

common/phitt/  ph (97) ,phd(97) , phb(97) ,phm(97)  0000080 

common/denst/  fd(97) ,fdd(97) , fdb(97) ,fdm(97)  0000090 

0000100 

output  the  hit  probabilities  O000110 

call  outfi (' fp.php' , ph,phd,phb,phm)  0000120 

0000130 

output  the  fragment  densities  0000140 

call  outfi (' fp.pfd' , fd, fdd, fdb, fdm)  0000150 

0000160 

return  0000170 

0000180 

end  0000190 

0000200 


****************************** ********************************  *********  O000210 
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0******** ************ ******************************************** ******* 
c 

subroutine  outfi (fname, v, vd, vb, vm) 
c 

c  outputs  hit  probabilities  or  fragment  densities 

c 

c  common  blocks 

common/repls/  nrepsr  fnreps 
common/donor/  nunits, htstft , htbaft 
common /rangs/  nrsegs, nrsegp 

common/ range/  rngmax,  rngseg,  rngsgh, rngscl, rngsix,  rngmin 
common/arang/  rangem ( 97 ) , ahf rn ( 97 ) , ahtop ( 97 ) 
common/names/  dname,  aname 
common/acptr/  htacft, wdacft, dpacft, htacfm 

local  arrays 

dimension  v(97) , vd(97) , vb (97) f  vm(97) 


c 

c 


c 

c 

c 

c 


type  declarations 

character  fname*6, dname*40, aname*40 

open (3, file=fname, status=' unknown' ) 

write (3, '  (3a) ' )  char (39) , dname, char (39) 

write (3, ' (2i3, Ip2el6 . 8) ' )  nrepsf nunits, htstft, htbaft 

write (3, ' (3a) ' )  char (39) , aname, char (39) 

write (3f' (Ip3el6.8)')  htacftr wdacftr dpacft 

write (3, ' (Ip2el6 . 8) ' )  rngmin, rngmax 

do  100  irange=l, nrsegs 

write (3, ' (Ip5el6 . 8) ' )  rangem (i range) , vd(irange) , 

+  vb(irange)  ,  vm(irange) ,v(irange) 

100  continue 

close (3) 

return 

end 


0000010 

0000020 

0000030 

0000040 

0000050 

0000060 

0000070 

0000080 

0000090 

0000100 

OQOOllO 

0000120 

0000130 

0000140 

0000150 

0000160 

0000170 

0000180 

0000190 

0000200 

0000210 

0000220 

0000230 

0000240 

0000250 

0000260 

0000270 

0000280 

0000290 

0000300 

0000310 

0000320 

0000330 

0000340 

0000350 

0000360 

0000370 

0000380 

0000390 


0***********************************************************************  0000400 
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q*********************************************************************** j^q  00010 


c 

R000020 

function  rnval (valmin, valmax, rndnxt) 

R000030 

c 

R000040 

c 

generates  a  random  value  between  two  specified  values 

R000050 

c 

R000060 

c 

type  declaration 

R000070 

double  precision  rndnxt 

R000080 

c 

R000090 

rnval=valmin+  (valmax- valmin)  *rndom (rndnxt) 

R000100 

c 

R000110 

return 

R000120 

c 

R000130 

end 

R000140 

c 

R000150 

c****** ********************************** ******************** ***********  R000160 

c** ************************************************************** *******  rooooiq 

c 

R000020 

function  rndom(uix) 

R000030 

c 

R000040 

c 

generates  random  numbers 

R000050 

c 

R000060 

c 

type  declarations 

R000070 

double  precision  ua, up, uix, uiy, ubl5, ub!6f 

R000080 

+  uxhi, uxalo,ulftlo, ufhi,  uk 

R000090 

c 

R000100 

c 

data  statement 

R000110 

data  ua/16807.0d0/,ubl5/32768.0d0/. 

R000120 

+  ubl6/65536. OdO/, up/2147483647. OdO/ 

R000130 

c 

R000140 

if (uix.eq. 0)  uix=uiy 

R000150 

100  uxhi=uix/ubl6 

R000160 

uxhi=uxhi-dmod  (uxhi,  1 .  OdO ) 

R000170 

uxalo= (uix-uxhi*ubl6) *ua 

R000180 

ulftlo=uxalo/ubl6 

R000190 

ulftlo=ulftlo-dmod (ulf tlof 1 . OdO ) 

R000200 

ufhi=uxhi*ua+ulftlo 

R000210 

uk=ufhi/ubl5 

R000220 

uk=uk-dmod (ukf 1 . OdO) 

R000230 

uix= ( ( (uxalo-ulftlo*ub!6) -up) + (ufhi-uk*ubl5) *ubl6) +uk 

R000240 

if (uix . It . 0 . 0 )  uix=uix+up 

R000250 

uiy=uix 

R000260 

rndom=int ( (uix*4 . 656612875d-10) *1.0d06) /1.0d06 

R000270 

if (rndom.eq. 0 . 0)  go  to  100 

R000280 

c 

R000290 

return 

R000300 

c 

R000310 

end 

R000320 

c  R000330 

c** ****************************************************** ***************  RQ00340 
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oooo  ooooooooo 


0***********************************************************************  vOO 00 10 


c  V000020 

function  vcr jr (abjr, c jr, hmm, fdmm)  V000030 

c  V000040 

c  computes  the  Jacobs-Roslund  critical  velocity  V000050 

c  V000060 

c  common  block  VO  00 070 

common/ cnstm/  qpi ,  hpi ,  pi ,  tpi ,  raddeg,  degrad,  smlrad,  smlnum,  bignum  V0 0 0  0 8 0 
c  V000090 

if (fdmm. eq. 0.0)  then  V000100 

vcrjr=bignum  V000110 

else  V000120 

vcr jr=abjr* (1 . 0+c jr*hmm/fdmm) / sqrt (fdmm)  V000130 

endif  V000140 

c  V000150 

return  V000160 

c  V000170 

end  V000180 

V000190 


***********************************************************************  y q 00200 


************************************** *********************************  V000010 


V000020 

function  vcrth (tc, ath,bth,hcm, facm2, fmgm)  V000030 

V000040 

computes  the  THOR  critical  velocity  V000050 

V000060 

common  block  V0 00 070 

common/ cnstm/  qpi, hpi, pi,  tpi,  raddeg, degrad,  smlrad,  smlnum, bignum  V0 00 080 

V000090 

if (fmgm.eq. 0 .0)  then  V000100 

vcrth=bignum  V0 00 110 

else  V000120 

vcrth=tc* ( (hcm*facm2) **ath) * (fmgm**bth)  V000130 

endif  V000140 

V000150 

return  VOO 0160 

V000170 

end  V000180 

V000190 


**********************************************  *************************  V000200 


q***********************************************************************  E000010 


c  E000020 

function  expon(arg)  E000030 

c  E000040 

c  computes  the  exponential  function  E000050 

c  limited  for  large  negative  arguments  E000060 

c  E000070 

c  common  block  E000080 

common/cnstm/  qpi, hpi, pi, tpi,  raddeg, degrad,  smlrad,  smlnum, bignum  E000090 
c  E000100 

expon=max (exp (arg) , smlnum)  E000110 

c  E000120 

return  E000130 

c  E000140 

end  E000150 

c  E000160 


c****************************** *****************************************  £000170 
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APPENDIX  C: 

LETHALITY  AND  VULNERABILITY  DATA 
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Intentionally  left  blank. 


78 


The  data  required  for  M107  and  TOW-2A  lethality  and  vulnerability  descriptions  is  summarized  in 
Tables  C-l  through  C-5. 

Table  C-l.  Weapon  Dimensions  and  Materials 


Component 

M107 

TOW-2A 

Weapon: 

Outside  Diameter 

6.10  in 

6.187  in 

Length 

23.90  in 

60.000  in* 

Warhead: 

Outside  Diameter 

6.10  in 

5.850  in 

Length 

23.90  in 

25.570  in 

Casing  Material 

steel 

aluminum 

Casing  Thickness 

0.63  in 

0.050  in 

Explosive  Type 

Comp.  B 

LX- 14 

Charge  Diameter 

4.84  in 

5.750  in 

Charge  Length 

20.30  in 

9.574  in 

Rocket  Motor. 

Outside  Diameter 

N/A 

5.837  in 

Length 

N/A 

18.700  in 

Casing  Material 

N/A 

steel 

Casing  Thickness 

N/A 

0.055  in 

Propellant  Type 

N/A 

GCV 

Charge  Diameter 

N/A 

5.782  in 

Charge  Length 

N/A 

6.320  in 

Container 

Material 

N/A 

steel 

Outside  Diameter 

N/A 

10.000  in 

Thickness 

N/A 

0.040  in 

Pallet: 

Orientation 

vertical 

horizontal 

Arrangement  (HxWxD) 

1x4x2 

3x1x4 

Weapon  Spacing 

1.20  in 

1.000  in* 

Elevation 

8.10  in 

8.100  in* 

Height  (with  pallet) 

32.00  in 

40.220  in 

Width 

28.00  in 

60.000  in 

Depth 

13.40  in 

43.160  in 

* 

estimated  value 
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Table  C-2.  Energetic  Material  Performance  Parameters 


Material 

Density 

(S/cm2) 

Detonation 

Pressure 

(kbar) 

Gurney 

Constant 

(m/s) 

Product 

Gases 

(mole/g) 

Molecular 

Weight 

Heat  of 
Detonation 
(cal/g) 

LX- 14 

1.83 

370.0 

2,948.0 

Not 

Required 

Not 

Required 

Not 

Required 

GCV  Propellant 

1.8 

Unavailable 

Unavailable 

0.042156 

24.635 

890.0 

Composition  B  3.065  0.0  1.70 

LX- 14* _ 3^6 _ 0.0  1.59 

GCV  Propellant*  4.8  0.0  1.3 


estimated  values 


Table  C-4.  THOR  Velocity  Equation  Constants 


Material 

n 

■91 

cv 

<K 

ev 

mild  steel 

3.6901 

0.889 

-0.945 

0.019 

1.262 

aluminum 

3.9356 

1.029 

-1.072 

-0.139 

1.251 
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APPENDIX  D: 

DONOR  FRAGMENTATION  GENERATOR  LISTING 
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INTENTIONALLY  LEFT  BLANK. 
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c************ ******************************************** *************** 
c 

program  dogen 
c 

C  DONOR  FRAGMENTATION  DATA  GENERATOR 

c  (Horizontal  Storage) 

c 

c  common  blocks 

common/frgwh/  iazwh(2500) , fmwhgm(2500) 
common/ f  rgrm/  iazrm(2500 )  ,  fmrmgm(2500) 
c 

c  type  declaration 

double  precision  rndnxt 
c 

c  define  constants 

call  const 
c 

c  read  and  compute  warhead  parameters 

call  elein (' wh.inp' , nfrwh, fmavwh, vfwhfs, plzwh, fmnwh) 
c 

c  read  and  compute  rocket  motor  parameters 

call  elein  ('  rm.  inp'  ,nfrrm,  fmavrm,  vf  rmf  s,plzrm,  fmnrm) 
c 

c  input  the  randomization  seed 

call  rndin (rndnxt) 
c 

c  generate  the  warhead  fragment  distriubtion 

call  f rdst  ( rndnxt ,  nfrwh,  fmavwh,  fmnwh,  nazwh,  If rwh,  iazwh,  fmwhgm, 

+  iazwhx) 

c 

c  generate  the  rocket  motor  fragment  distriubtion 

call  frdst  ( rndnxt ,nfrrm,  fmavrm,  fmnrm,  nazrm,  lfrrm,  iazrm,  fmrmgm, 

+  iazrmx) 

c 

c  generate  fragmentation  data  file 

call  mkf il (plzwh, nazwh, If rwh, nfrwh, vfwhfs, fmnwh, iazwhx, 

+  plzrm,  nazrm,  lfrrm,  nf rrm,  vfrmfs,  fmnrm,  iazrmx) 

c 

stop  '  9 
c 

end 

c 

c** ********************************************************** *********** 


D000010 

D000020 

D000030 

D000040 

D000050 

D000060 

D000070 

D000080 

D000090 

D000100 

D000110 

D000120 

D000130 

D000140 

D000150 

D000160 

D000170 

D000180 

D000190 

D000200 

D000210 

D000220 

D000230 

D000240 

D000250 

D000260 

D000270 

D000280 

D000290 

D000300 

D000310 

D000320 

D000330 

D000340 

D000350 

D000360 

D000370 

D000380 

D000390 

D000400 

D000410 

D000420 

D000430 


c** **************** ******** ********************************** ***********  COO 0010 


C  C000020 

subroutine  const  C000030 

c  C000040 

c  defines  constants  C000050 

c  C000060 

c  common  block  C000070 

common/ cnstm /  pi, twth,  dazdeg, idzdeg, azxdeg  CO 0 0 0 8 0 

c  C000090 

pi=4 . 0*atan (1 . 0)  C000100 

twth=2 . 0/3 . 0  C000110 

dazdeg=10 .0  CO 00 120 

idzdeg=int (dazdeg)  C000130 

azxdeg=110 • 0  C000140 

c  C000150 

return  CO 00 160 

c  C000170 

end  C000180 

c  C000190 


C***************** *********************************************** *******QQQQ  2Q0 
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oo  oo  oo  oo  oo  oo  oo  oo  oo  oo  oo  oo 


c********** ****************** *******************************************  E000010 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  elein  (inf ile, nfrag,  frmav, vf f s, plzdeg, f nrnnin) 
reads  and  computes  cylindrical  element  parameters 
common  block 

common/ cnstm/  pi, twth, dazdeg, idzdeg, azxdeg 

type  declaration 
character*  6  infile 

open  data  file 

open (3, f ile=inf ile, status='  old' ) 
read  casing  data 

read  (3,  *)  rhea,  romm,  drmm,  htmm,plzdeg,  f  rmmin 

read  fill  data 
read(3, *)  rhfl, pdfl, vgfl 
if (pdfl.eq.0.0)  read (3,*)  fln,flm,flq 

close  data  file 
close (3) 

convert  dimensions  to  centimeters 
rocm=0 . l*romm 
drcm=0 .  l*drmm 
htcm=0.1*htmm 

compute  the  casing  thickness  and  inside  diameter  in  inches 
ar  in=drmm/  2  5 . 4 
diin=  ( romm-drmm)  /25 . 4 

compute  the  inside  radius  in  centimeters 
r  i  cm=  r  o  cm-dr  cm 
ricm2=ricm*ricm 

compute  the  casing  mass  in  grams 
camgm=rhca*pi*  (rocm*rocm-ricm2)  *htcm 
write  (*, '  (a, lpe!2 .5, a) ' )  ' 


casing  mass  =',camgm. 


compute  the  fill  mass  in  grams 
flmgm=rhfl*pi*ricm2*htcm 
write (*,' (a,lpel2.5,a) 9  )  9 


gf 


compute  the  fragment  velocity 
vfms^vgf  1/sgrt  (camf  1+0 .5) 
vffs=3 .2808 J99*vfms 
write (*, 9  (a, lpel2 .5, a) ' )  ' 


fragment  velocity  =',vfms,'  m/s' 


compute  the  average  fragment  mass  in  grams 
f  rmav— 


E000020 

E000030 

E000040 

E000050 

E000060 

E000070 

E000080 

E000090 

E000100 

E000110 

E000120 

E000130 

E000140 

E000150 

E000160 

E000170 

E000180 

E000190 

E000200 

E000210 

E000220 

E000230 

E000240 

E000250 

E000260 

E000270 

E000280 

E000290 

E000300 

E000310 

E000320 

E000330 

E000340 

E000350 

E000360 

E000370 

E000380 

E000390 

E000400 

E000410 

E000420 

E000430 

E000440 

E000450 

E000460 

E000470 

E000480 

E000490 

E000500 

E000510 

E000520 

E000530 

E000540 


fill  mass  =',flmgm, 

compute  the,  ratio  of  casing  to  fill  mass 
camfT=camgm/  f  lmgm 

compute  the  detonation  pressure  if  necessary 
if  •?<?•  °4  °>  1=15  *  58*rhf  l*rhf  l*f  ln*sqrt  (f  lm*f  lq) 

write (*,' (a, lpel2 .5, a)  )  '  detonation  pressure  =',pdfl,'  kbar' 

compute  the  Gurney  constant  if  necessary 

if  j[vgfl.eg.  0.0)  vgf 1=233 . 0*sqrt  (pdfl)  /rhf  1**0 . 6  E000550 

write (*, 9  (a, lpel2 . 5, a) ' )  '  Gurney  constant  =', vgfl, '  m/s'  E000560 

E000570 
E000580 
E000590 
E000600 
E000610 
E000620 
E000630 

+  6 . 7 62 e2  * sort ( 1 . 0+0 . 5  *  camf  1 ) * (drin* ( (diin+drin) **1.5) /diin) /pdfl  E000  650 
wnte(*, '  (a,  lpel2.5,a) ' )  '  average  fragment  mass  =',frmav,'  g'  §000660 

size  if  needed) 


compute  the  number  of  fragments  (adjust  sector 
100  nfrag=int (camgm*azxdeg/ (360 . 0*frmav) ) 
if (nfrag.gt .2500)  then 

azxdeg=2500 . 0*azxdeg/f loat (nf rag) 
go  to  100 
enaif 


E000670 

E000680 

E000690 

E000700 

E000710 

E000720 

E000730 


write  *,  3, IS)  )  '  number  of  fragments  ='nf rag  E000740 

write (*,  (a, lpel2 .5, a/)  )  '  azimuthal  sector  size  =', azxdeg,'  deg'  E000750 

„  E000760 

return  E000770 

,  E000780 

end  E000790 

c************** ********************************************** ***********  E000810 


C 

c 
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o  o 


q**** ********************************************************  ***********  F000010 


c 

subroutine  frdst  (rndnxt, nf rag,  frmav,  frmmin,naz,  If rag,  iaz,  f  rmgm, 

+  iazmax) 

c 

c  generates  a  fragment  distriubtion 

c 

c  common  block 

common/cnstm/  pi, twth, dazdeg, idzdeg, azxdeg 

type  declaration 
double  precision  rndnxt 
c 

c  local  arrays 

dimension  iaz (2500) ,  frmgm(2500) 
c 

c  generate  fragment  masses  with  Mott  distribution 

c  and  azimutal  angles  with  even  distribution 

hf rmav=0 . 5*f rmav 
iazmax=0 
If rag=0 

do  100  ifrag=l,nfrag 
rnddum=  rndom  ( rndnxt ) 

f  rmgm  (if  rag)  =-hf  rmav*alog  ( rnddum*  rnddum) 

iaz (if rag) =idzdeg*int (azxdeg* rndom (rndnxt) /dazdeg) +idzdeg 
iazmax=max (iazmax, iaz (if rag) ) 
if  (f  rmgm  (if  rag)  .gt.frmmin)  If  rag=lf  rag+1 
100  continue 

naz=iazmax/ idzdeg 
c 

return 

c 

end 

c 

Q**  ****************  **********************  ************************  *★★*★*•* 


F000020 

F000030 

F000040 

F000050 

F000060 

F000070 

F000080 

F000090 

F000100 

F000110 

F000120 

F000130 

F000140 

F000150 

F000160 

F000170 

F000180 

F000190 

F000200 

F000210 

F000220 

F000230 

F000240 

F000250 

F000260 

F000270 

F000280 

F000290 

F000300 

F000310 

F000320 

F000330 

F000340 

F000350 
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O  O  O  O  oooooooooooo  ooo  oooo  oo 


***********************************************************************  M000010 


subroutine  mkf  il  (plzwh,  nazwh.  If rwh,  nf  rwh,  vfwhfs,  fmnwh,  iazwhx, 
+  plzrm,  nazrm,  if  rrm,  nf  rrm,  vf  rmf  s,  fmnrm,  iazrmx) 

generates  fragmentation  data  file 

common  blocks 

common/ cnstm/  pi, twth, dazdeg, idzdeg, azxdeg 
common /frgwh/  iazwh(2500) , fmwhgm(2500) 
common/f rgrm/  iazrm(2500) , fmrmgm(2500) 

output  fragmentation  data  with 
azimuthal  and  polar  angles  interchanged 
naz=max (nazwh, nazrm) 
iazmax=max (iazwhx, iazrmx) 

If rag=lf rwh+lf rrm 

open (3, f ile=' f r .out' , status=' unknown' ) 
write (3,' (3x, Ip2el5 .7) ' )  dazdeg, max (plzwh, plzrm) 
write (3, ' (3x, 3i5) ' )  naz,lfrag, 1 
write (4, lpel5.7)  0.04 

write (4, Ip4el5 .7)  3.6901,0.889,-0.945,0.019 
write (4, Ip4el5 .7)  -2.478,0.138,0.835,0.761 
write(3, ' (3x,lp4el5.7) ' )  (0 . 05, iz=l, naz) 
do  100  iwrite=idzdeg, iazmax, idzdeg 

call  adf rg (iwrite, nf rwh, fmwhgm, iazwh, fmnwh, vfwhfs) 
call  adf  rg  (iwrite,  nf  rrm,  fmrmgm,  iazrm,  fmnrm,  vf  rmfs) 

100  continue 
close (3) 

return 

end 


M000020 
M000030 
M000040 
M000050 
M000060 
M0  00 070 
M000080 
M000090 
M000100 
M000110 
M000120 
M000130 
M000140 
M000150 
M000160 
M000170 
M000180 
M000190 
M000200 
M000210 
M000220 
M000230 
M000240 
M000250 
M000260 
M000270 
M000280 
M000290 
M000300 
M000310 
M000320 
M000330 
M000340 


**************************************************************  *********  MOO  0350 


***********************************************************************  a00 0010 


subroutine  adfrg  (iwrite,  nf rag,  frmgm,  iaz,  fmn,  vffs) 


common  blocks 

common/ cnstm/  pi, twth, dazdeg, idzdeg, azxdeg 
local  arrays 

dimension  iaz (2500) , frmgm(2500) 

do  100  if rag=l,nfrag 

if (f rmgm(if rag) .gt.fmn)  then 
if (iaz (if rag) .eq. iwrite)  then 
acm2=0 . 5199*f rmgm(if rag) **twth 
ain2-0 . 1550003*acm2 
frmlb=2 .2046226e-3*f rmgm(if rag) 
frmgr=7000 .0*frmlb 
amin lb=a in2 / f rmlb 

write  (3, ' (i3, Ip4el5 . 7) ' )  iaz (if rag) , frmgr,  aminlb, vffs, 1 . 5 
endif 
endif 

100  continue 
return 
end 


A000020 
A000030 
A000040 
A000050 
A000060 
A000070 
A0  00 080 
A000090 
A000100 
A000110 
A000120 
A000130 
A000140 
A000150 
A000160 
A0  00 170 
A000180 
A000190 
A000200 
A000210 
A000220 
A000230 
A000240 
A000250 
A000260 
A000270 
A000280 


***********************************************************************  A00 0290 
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c************** ************************** *******************************  R000010 


c 

R000020 

subrout ine  rndin ( rndnxt ) 

R000030 

c 

R000040 

c 

obtains  randomization  seed  input 

R000050 

c 

R000060 

c 

type  declarations 

R000070 

double  precision  rndinp, rndnxt 

R000080 

c 

R000090 

100  write (*, f  (a\) f )  f  Enter  Monte  Carlo  Seed  (1  to  2147483646):  9 

R000100 

read(*, *,err=100)  rndinp 

R000110 

if (rndinp. It. 1.0. or .rndinp. gt . 2147483646.0)  go  to  100 

R000120 

c 

R000130 

rnddum=rndom  ( rndinp) 

R000140 

rndnxt =0 . OdO 

R000150 

c 

R000160 

return 

R000170 

c 

R000180 

end 

R000190 

c 

R000200 

c** ****************************************************** ***************  RQ00210 


0** ****** ****************** ******** ****************** ****** ************* 

c 

function  rndom(uix) 
c 

c  generates  random  numbers 

c 

c  type  declarations 

double  precision  ua, up, uix,uiy, ub!5, ubl6, uxhi, uxalo, ulftlo, ufhi, 

+  uk 

c 

c  data  statement 

data  ua/16807.0d0/,ubl5/32768.0d0/,ubl6/ 65536. 0d0/f 
+  up/2147483647 . OdO/ 

c 

if ( uix . eq . 0 )  uix=uiy 
100  uxhi=uix/ubl6 

uxhi=uxhi~dmod (uxhi, 1 . OdO) 

uxalo= (uix-uxhi*ubl6) *ua 

ulftlo=uxalo/ubl6 

ulftlo^ulftlo-dmod (ulftlo, 1 . OdO) 

ufhi=uxhi*ua+ulftlo 

uk=ufhi/ubl5 

uk=uk-dmod (uk, 1 . OdO ) 

uix=( ( (uxalo-ulftlo*ubl6)-up)+(ufhi-uk*ubl5) *ubl6)+uk 

if (uix . It . 0 . 0 )  uix=uix+up 

uiy=uix 

rndom=int ( (uix*4 . 656612875d-10) *1 . 0d06) /I . 0d06 
if (rndom.eq. 0 . 0)  go  to  100 
c 

return 

c 

end 

c 

c** ******************************************** **********  *************** 


R000010 
R000020 
R000030 
R000040 
R000050 
R000060 
R000070 
R000080 
R000090 
R000100 
R000110 
R000120 
R000130 
R000140 
R000150 
R000160 
R000170 
R0 00180 
R000190 
R000200 
R000210 
R000220 
R000230 
R000240 
R000250 
R000260 
R000270 
R000280 
R000290 
R000300 
R000310 
R000320 
R000330 
R000340 
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Intentionally  left  blank. 


8S 


NO.  OF 

COPIES  ORGANIZATION 


2  ADMINISTRATOR 

DEFENSE  TECHNICAL  INFO  CTR 
ATTN  DTIC  DDA 
CAMERON  STATION 
ALEXANDRIA  VA  22304-6145 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
ATTN  AMSRL  OP  SD  TA 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1145 

3  DIRECTOR 

US  ARMY  RESEARCH  LAB 
ATTN  AMSRL  OP  SD  TL 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1145 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
ATTN  AMSRL  OP  SD  TP 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1145 


ABERDEEN  PROVING  GROUND 

5  DIR  USARL 

ATTN  AMSRL  OP  AP  L  (305) 


NO.  OF 

COPIES  ORGANIZATION 


NO.  OF 

COPIES  ORGANIZATION 


1  DIRECTOR 

USA  TECH  CTR  FOR  EXPLOSIVE  SAFETY 
ATTN  SMCAC  ESL  C  DOYLE 
SAVANNA  IL  61074-9639 

1  DIRECTOR 

USA  TECH  CTR  FOR  EXPLOSIVE  SAFETY 
ATTN  SMCAC  DO  N  MCCULLOUGH 
SAVANNA  IL  61074-9639 

3  US  ARMY  WATERWAYS  EXP  STATION 
ATTN  CEWES  SE 

K  DAVIS 
M  FORD 
C  JOACHIM 
3909  HALLS  FERRY  RD 
VICKSBURG  MS  39180-6199 

2  CHAIRMAN 

DOD  EXPLOSIVE  SAFETY  BOARD 
ATTN  DDESB  KT 
J  WARD 
C  CANADA 

2461  EISENHOWER  AVE 
ALEXANDRIA  VA  22331-0600 

4  DEFENSE  AMMUNITION  LOGISTICS  ACTIVITY 
ATTN  AMCPM  AL 

R  ROSSI 
RHO 
E  GOON 
B  WILLIAMSON 

PICATINNY  ARSENAL  NJ  07806-5000 

1  COMMANDER 

ATTN  SMCAR  AEE  W 
PLU 

US  ARMY  ARDEC 

PICATINNY  ARSENAL  NJ  07806-5000 

1  COMMANDER 

ATTN  SMCAR  AEM 
B  WILLIAMSON 
US  ARMY  ARDEC 

PICATINNY  ARSENAL  NJ  07806-5000 

1  COMMANDER 

ATTN  SMCAR  ASA 
D  MILLER 
US  ARMY  ARDEC 

PICATINNY  ARSENAL  NJ  07806-5000 


1  OFFICER  IN  CHARGE 

WHITE  OAK  LABORATORY  NSWC 
ATTN  R15  M  SWISDAK 
10901  NEW  HAMPSHIRE  AVE 
SILVER  SPRING  MD  20902-5000 

2  NAVAL  FACILITIES  ENGINEERING  SERV  CTR 
ATTN  CODE  LSI 

J  TANCRETO 
RMURTHE 

PORT  HUENEME  CA  93043 

1  DIRECTOR 

EXPLOSIVES  HAZARD  REDUCTION  DIR 

ATTN  ASC  YOCO  EHR 

JJENUS 

EGLIN  AFB  FL  32542 

2  US  ARMY  CORPS  OF  ENGINEERS 
ATTN  CEHND  ED  CS 

R  WRIGHT 
RHASSE 
PO  BOX  1600 

HUNTSVILLE  AL  35807-4301 

ABERDEEN  PROVING  GROUND.  MD 

13  DIR,  USARL 

ATTN:  AMSRL-WT-T,  W.  MORRISON 
AMSRL-WT-TB, 

P.  BAKER 
R.  FREY 
O.  LYMAN 

J.  WATSON 
F.  GREGORY 

V.  BOYLE 

W.  HILLSTROM 
E.  MCDOUGAL 
W.  LAWRENCE 

K.  BENJAMIN 
T.  DORSEY 

AMSRL-SL-BV,  J.  PLOSONKA 

2  DIR,  USA  AMSAA 
ATTN:  AMXSY-GC, 

E.  CHRISTMAN 
A.  WONG 


90 


USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  publishes.  Your  comments/answers 
to  the  items/questions  below  will  aid  us  in  our  efforts. 

1.  ARL  Report  Number  ARL-TR-949 _ Date  of  Report  January  1996 _ 

2.  Date  Report  Received _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or  other  area  of  interest  for  which  the  report 

will  be  used.) _ 


4.  Specifically,  how  is  the  report  being  used?  (Information  source,  design  data,  procedure,  source  of  ideas,  etc.) 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or  dollars  saved,  operating  costs 
avoided,  or  efficiencies  achieved,  etc?  If  so,  please  elaborate. 

'  —— 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate  changes  to 
organization,  technical  content,  format,  etc.)  _ _ _ 


Organization 


CURRENT  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  Current  or  Correct  address  above  and  the 
Old  or  Incorrect  address  below. 


Organization 


OLD  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 

(Remove  this  sheet,  fold  as  indicated,  tape  closed,  and  mail.) 
(DO  NOT  STAPLE) 
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FIRST  CLASS  PERMIT  NO  0001,APG,MD 


POSTAGE  WILL  BE  PAID  BY  ADDRESSEE 


DIRECTOR 

U.S.  ARMY  RESEARCH  LABORATORY 
ATTN:  AMSRL-WT-TB 

ABERDEEN  PROVING  GROUND,  MD  21005-5066 


NO  POSTAGE 
NECESSARY 
IF  MAILED 
IN  THE 

UNITED  STATES 


